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Abstract 


The present thesis addresses the design of structure-preserving numerical methods that 
emanate from the general equation for non-equilibrium reversible-irreversible coupling 
(GENERIC) formalism. 


For conservative mechanical systems, the GENERIC reduces to the purely Hamiltonian 
case. Considering this special case first, novel energy-momentum (EM) consistent time- 
stepping schemes in the realm of molecular dynamics are introduced that take periodic 
boundary conditions, three-body potentials and interatomic functional potentials into 
account. This method is thermodynamically consistent because it preserves the energy 
and the symmetries of the system, independent of the time-step size. 


Afterwards, thermodynamical systems are considered in the framework of the 
GENERIC, and a new variational formulation for large-strain thermoelasticity is pro- 
posed. This variational formulation emanates from the GENERIC framework and allows 
for the free choice of the thermodynamic state variable among three options: (i) the 
absolute temperature, (ii) the internal energy density, or (iii) the entropy density. In this 
context the notion "GENERIC-consistent space discretization" is introduced, facilitating 
the design of novel energy-momentum-entropy (EME) consistent schemes. This method 
is thermodynamically consistent in the sense that the solution strictly obeys the first and 
second law of thermodynamics and the symmetries of such a system, independent of the 
size of the time-step. 


Finally, large-strain thermo-viscoelasticity is considered in the context of the 
GENERIC framework. A mixed finite element approach for the discretization in space is 
proposed that incorporates a GENERIC-consistent space discretization. Depending on 
the choice of the thermodynamic state variable, the plain mid-point rule already yields 
partially structure-preserving schemes on its own. Thus, this newly developed GENERIC- 
based weak form is particularly well suited for the design of structure-preserving methods. 
In all cases, numerical investigations are presented that confirm the theoretical findings 
and shed light on the numerical stability of the newly developed schemes. 


Keywords: structure-preserving space-time integration, GENERIC, molecular dynamics, 
periodic boundary conditions, interatomic potentials, thermoelasticity, thermo-viscoelas- 
ticity, finite element method, finite deformation 


Zusammenfassung 


Die vorliegende Dissertation behandelt die Entwicklung von strukturerhaltenden nu- 
merischen Methoden, welche aus dem GENERIC (Akkronym für general equation for 
non-equilibrium reversible-irreversible coupling) Formalismus hervorgehen. 


Für konservative mechanische Systeme reduziert sich das GENERIC auf die hamiltoni- 
sche Beschreibung. Zunächst werden für diesen Spezialfall neuartige energy-momentum 
(EM) konsistente Zeitschrittverfahren für den Bereich molekulardynamischer Systeme 
eingeführt, welche periodische Randbedingungen, Dreikórperpotentiale und interatomare 
Funktionalpotentiale berücksichtigen. Eine derartige Methode ist thermodynamisch kon- 
sistent, da sie die Energie und die Symmetrien einer solchen Struktur unabhàngig von der 
Zeitschrittgröße bewahrt. 


Anschließend wird der Fall von thermodynamischen Systemen im GENERIC Formalis- 
mus betrachtet und eine neue Variationsformulierung für die Thermoelastizität unter 
Berücksichtigung großer Deformationen vorgeschlagen. Diese Variationsformulierung 
entstammt dem GENERIC-Gerüst und ermóglicht die freie Wahl der thermodynamischen 
Zustandsvariable unter folgenden drei Optionen: (i) der absoluten Temperatur, (ii) der 
inneren Energiedichte, oder (iii) der Entropiedichte. In diesem Zusammenhang wird der 
Begriff der “GENERIC-konsistenten ráumlichen Diskretisierung" eingeführt, welcher den 
Entwurf neuartiger energy-momentum-entropy (EME) konsistenter Schemata ermóglicht. 
Eine solche vorgeschlagene Methode ist thermodynamisch konsistent, da die Lósung 
streng dem ersten und zweiten Hauptsatz der Thermodynamik und den Symmetrien eines 
solchen Systems folgt, unabhängig von der Größe des Zeitschritts. 


Schließlich wird die Thermo-Viskoelastizität unter Berücksichtigung großer Deforma- 
tionen im Kontext des GENERIC-Gerüsts betrachtet. Für die Diskretisierung im Raum 
wird ein Ansatz mit gemischten finiten Elementen vorgeschlagen, welcher sich einer 
GENERIC-konsistenten ráumlichen Diskretisierung bedient. Abhángig von der Wahl der 
thermodynamischen Variable liefert die einfache Mittelpunktregel bereits teilweise struk- 
turerhaltende Schemata. Damit eignet sich die neuartige GENERIC-basierte schwache 
Form ideal für die Entwicklung strukturerhaltender Verfahren. In allen Fállen werden 
numerische Untersuchungen vorgestellt, die die theoretischen Erkenntnisse bestätigen 
und Aufschluss über die numerische Stabilität der neu entwickelten Verfahren geben. 


ili 


Zusammenfassung 


Schlüsselwörter: Strukturerhaltende Raum-Zeit Integration, GENERIC, Molekulardy- 
namik, periodische Randbedingungen, interatomare Potentiale, Thermoelastizitát, Ther- 
mo-Viskoelastizitat, Finite Elemente Methode, große Deformationen 
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1. Introduction! 


Mathematical models are extensively used by engineers to gain insight into the quantitative 
and structural behavior of a system. Although simple systems can be understood and 
analyzed more easily, complexity arises when a more accurate representation of reality is 
required, necessitating the replacement of pencil and paper by large computer simulations 
in favor of more powerful design tools. 


At the heart of every dynamical computer simulation lies the numerical method, also 
known as the "integrator", which constructs an approximation of the exact solution by 
providing discrete snapshots of the system over time for a given problem. Today, driven 
by the exponential growth of computational power, the qualitative properties of the inte- 
grator itself become vital to the success of the simulation, as longer simulation durations 
require more stable numerical methods. While the equations of motion for a dynamic 
system may "hide" important properties of the physical system, such as conservation laws, 
these properties are not inherited by the numerical method in general, which can cause 
unphysical time evolution and thus introduces a source for numerical instabilities. One 
feasible approach that counters this physical inconsistency is the use of subtle numer- 
ical methods that are numerically accurate and capable of reliably capturing the main 
characteristics of the underlying physical system. This motivated the development of 
structure-preserving numerical methods, also known as “geometric integrators”, for a 
wide array of fields in applied and theoretical science. Geometric integrators are designed 
to respect the fundamental physics of the underlying mathematical model, ranging from 
time reversibility to first integrals, by preserving the geometric properties in a discrete 
setting. Most attention has been focused hereby on Hamiltonian or Lagrangian systems 
due to their well-understood mathematical (geometric) structure. Based on the common 
underlying structure, structure-preserving methods can be easily generalized from simple 
systems such as spring-mass arrays to complex nonlinear elastodynamical structures. 
Therefore, reports of excellent long-term behavior and numerical stability for a huge 
variety of systems have laid the foundation for the success of this method. 


1 This chapter is based on the introductions given in [1-4] 
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11. Structure-preserving integration 


Structure-preserving numerical methods have played a major role in the development of 
numerical time integration methods for decades, and monographs such as [5-7] beautifully 
demonstrate the advancement of this field. While numerical integrators have traditionally 
been viewed solely as an approximation to a time-continuous problem, geometric integra- 
tors consider the integration scheme itself as a discrete dynamical system with its own 
characteristics, including balance laws, symmetries, symplecticity, reversibility, and many 
more. Even though it was developed within a different context, the well-known mid-point 
rule can be viewed as a member of this general class because it preserves the symmetry 
and the symplecticity. 


Different classes of structure-preserving numerical methods can be distinguished by the 
features inheritated from the continuous problem and the applied numerical techniques. 
The earliest developments of structure-preserving numerical methods can be traced back 
to symplectic methods, first introduced in [8]. Early schemes of this kind were explicit 
symplectic integrators [9] and have been applied to various problems of celestial mechanics 
(see, e.g., [10, 11]), followed by implicit schemes [12, 13] based on [14]. Later, a symplectic 
family of Runge-Kutta schemes were constructed independently in [15, 16]. 


For mechanical systems, most researchers have focused their attention on momenta, 
energy, and symplectic structures. Unfortunately, simultaneously preserving the mo- 
menta, energy and symplecticity for a fixed time-step method is not possible (for more 
details, see [17]). Therefore, many geometric integrators in the realm of Lagrangian and 
Hamiltonian mechanics can be gathered around two classes: variational integrators and 
energy-momentum (EM) integrators. 


1.1.1. Variational integrators 


A more contemporary approach towards symplectic methods is based on a variational 
nature and is thus termed *variational integrator". Variational integrators preserve the 
symplecticity and the momenta of the underlying system and, due to its variational nature, 
allows the extension to non-conservative systems by utilizing the discrete counterpart of 
the Lagrange-d'Alembert principle. The origin of this method can be traced back to [18] 
and [19, 20]. Based on those concepts, a theory of discrete Lagrangian and Hamiltonian 
mechanics was provided in [21], which has been applied to mechanical problems with 
multi-symplectic geometry [22, 23], dynamical systems evolving on nonlinear manifolds 
[24, 25], structural elements [26-28], contact and impact problems [29], multibody dy- 
namics and control [30, 31], stochastic differential equations [32], constrained and forced 
problems [33] and dynamic of fluids [34] (see, e.g., [35] for more examples). 


Sharp phase portraits and long numerical simulations demonstrate the strong performance 
of this method, especially in explicit calculations. Attempts have been made to express 
thermoelasticity in the framework of Lagrangian/Hamiltonian mechanics by introducing 


1.1. Structure-preserving integration 


the concept of thermal displacement [36]. Unfortunately, expressing Fourier's law of 
heat conduction within this framework remains an unsolved issue [35]?. In addition, the 
loss of symplecticity and thus the loss of one key feature of variational integrators is a 
consequence of dissipation. For that matter, it seems that structure-preserving numerical 
methods for thermodynamical systems fit better within the energy-momentum integrator 
approach, which has been preferentially chosen for this work. Nevertheless, besides the 
preservation of the symplectic structure and the momenta, variational integrators are 
known for their improved long-term energy behavior? compared to classical methods and 
are excellent numerical time-stepping schemes. 


1.12. Energy-momentum integrators 


EM algorithms are an additonal, frequently used class of structure-preserving integration 
schemes which, in contrast to variational integrators (symplectic-momentum integrators), 
preserve the momenta and the energy of the system in the discrete setting. This method 
has its roots in [40-44], whereas the first energy-momentum scheme for nonlinear elasto- 
dynamics was proposed in [45]. Based on these previous reports, a systematic approach 
using the discrete gradient operator was developed in [46]. This second-order accurate 
operator ensures the preservation of the energy and the momenta by design and, due to 
its systemized construction, its remarkable robustness, and its good qualitative accuracy, 
has made the energy-momentum integrator a popular choice for simulating the governing 
equations of particle dynamics [40, 41], rigid bodies [43], nonlinear solid mechanics [45, 
47, 48], nonlinear shells [49-51] and rods [52-54] multi-body dynamics [55, 56], gradient 
systems [57], and general PDEs [58]. For a comprehensive overview of energy-momentum 
integrators in the context of Hamiltonian dynamics, see [59]. 


1.13. Dissipation 


The conservation of energy is not always a desired feature of Hamiltonian or Lagrangian 
systems because occasionally there is a need for dissipation in resolving high-frequency 
modes, which is often related to the spatial discretization of the problem. These "spurious" 
high frequencies are not present in real systems due to physical damping. However, physi- 
cal damping is seldom addressed in the mathematical modeling and is usually compensated 
for by introducing controllable numerical dissipation schemes. The first well-known rep- 
resentatives of such schemes can be found in the family of the Newmark method [60], the 


Efforts have been made in [37] to extend the Hamiltonian to account for the second sound phenomenon 
by adding a dissipative entropy-flux perturbation term and utilizing the Lagrange-D'Alembert principle. A 
similar approach has been followed for Fourier-type heat conduction in [38, 39]. It is not clear to the author if 
and how this adiabatic method can be generalized for arbitrary thermodynamical processes. 

For conservative systems, variational integrators capture the evolution of a “perturbated” Hamiltonian exactly. 
This modified energy-level set is, most likely, closely related to the true energy-level set and explains the 
characteristic energy oscillation of this method. 
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HHT method developed in [61], or the Wilson-0 method [62] among other methods that 
have been successfully implemented in common software packages. Unfortunately, these 
classical dissipation schemes, being developed in the linear regime, lose their dissipative 
character in the nonlinear regime [63]. In that sense, several energy-decaying numerical 
methods have been proposed [64-67]. 


Nevertheless, the philosophy of the current work is to incorporate physical dissipative 
processes into the mathematical model instead of considering numerical dissipation 
schemes. Most real structures can only be idealized as Hamiltonian or Lagrangian systems, 
yet are of relevance in many areas of science and engineering. 


Polymeric solids and rubber-like materials are a perfect example of structures that exhibit 
viscoelastic material features and are highly temperature sensitive, thus being subjected 
to viscous dissipation and dissipation due to the conduction of heat. Consequently, a 
large class of problems involving dissipative mechanisms can not be expressed within 
the framework of Hamiltonian or Lagrangian mechanics, including, for example, the 
aforementioned Fourier-type heat conduction, which arises in a wide range of industrial 
applications, including the following: 


* Civil engineering: Tuned mass damper (bridges, high-rise buildings), dissipative 
connections in high-rise buildings for seismic protection, and welding 


e Mechanical engineering: Shock absorbers (automotive and railway), brake discs, 
welding, and energy dissipation in vehicle crashes 


Unfortunately, structure-preserving integrators for thermodynamical systems have re- 
ceived little attention, one reason being the lack* of a well-understood geometric structure 
guiding the construction of such numerical method. Numerous attempts have been made 
to extend structure-preserving time-stepping schemes to the domain of non-conservative 
mechanical systems, such as port-Hamiltonian systems [68], viscoelasticity [69, 70], elasto- 
plasticity [71], thermoelasticity [72, 73], thermo-viscoelasticity [74-76], isothermal and 
non-isothermal fiber-reinforced continuas [77-80]. Unlike the Lagrangian/Hamiltonian 
case, these structure-preserving schemes do not emerge from a unifying theory. 


However, the metriplectic structure of the general equation for non-equilibrium reversible- 
irreversible coupling (GENERIC) formalism provides a double-generator framework that 
expresses thermomechanical models of dissipative materials (or generalized standard 
materials) in a unifying formalism. Its formulation is based on an additive decomposition 
of the time-evolution equations into a reversible part and a dissipative part. It can be seen 
as a natural extension of Hamiltonian mechanics to the dissipative regime (see, e.g., [81] 
for a relation to port-Hamiltonian systems). 


Originally developed in the context of complex fluids [82, 83], GENERIC has been applied 
to a vast range of problems (predominantly in the domain of complex fluids), such as 
the reptation model for entangled linear polymers [84], polymer blends [85], colloidal 


^ Atleast to the computational mechanics' community 
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suspensions [86], two-phase systems [87], relativistic hydromechanics [88], discrete for- 
mulations of hydromechanics for simulations [89], thermodynamically guided simulations 
[89], diffusion through polymeric membranes [90], rheological model of suspension of red 
blood cells [91], structure-preserving neural networks [92], and the theory of quantum 
systems [93], to mention only a few examples. We refer to the book by Ottinger [94] for a 
comprehensive account of the GENERIC formalism up to the year 2005 and to [95] for a 
huge list of advanced developments up to the year 2017. 


Focusing on thermodynamic models for non-isothermal solids, as one of the main themes 
of the present work, an early application of the GENERIC formalism to finite-strain thermo- 
elasticity is credited to [96], albeit in an Eulerian setting which is quite uncommon in solid 
mechanics, and to [97] for continuum damage mechanics. Later, the Lagrangian setting 
was used in [98, 99] and [100] to develop the GENERIC framework for non-isothermal 
solid mechanics. Moreover, [98, 99] preferred to use the absolute temperature as the 
thermodynamic state variable, while in [100] a special form of GENERIC is devised that 
makes the free choice of the thermodynamic state variable possible. 


In the field of computational solid mechanics, Romero [101, 102] recognized at an early 
stage the great potential of the GENERIC framework for the design of structure-preserving 
time-stepping schemes. Since the GENERIC framework automatically ensures the thermo- 
dynamic admissibility of the time-evolution equations, it provides an ideal starting point 
for the development of thermodynamically consistent (TC) integrators. The solutions of 
this integrator could be classified as thermodynamically guided simulations, in analogy 
to the thermodynamically guided simulation of consistent coarse-graining schemes (e.g., 
[103, 104]), as TC integrators comply with the first and the second law of thermodynamics, 
independent of the size of the time-step. Therefore, TC integrators may also be termed 
“energy-entropy” integrators. 


If TC integrators also respect the symmetries of the underlying mechanical system, they 
can be viewed as an extension to the dissipative regime of Energy-Momentum (EM) 
integrators previously developed for Hamiltonian systems with symmetry. The GENERIC 
framework also facilitates a concise characterization of momentum maps and associated 
conservation properties (see [105]). 


Hence, the GENERIC formalism provides a common structure for non-conservative sys- 
tems, similarly to Hamiltonian dynamics for conservative problems, and can therefore be 
used in a similar manner to guide the design of structure-preserving numerical schemes 
for thermodynamical problems. Consequently, GENERIC provides a solid theoretical 
foundation for the design of energy-momentum-entropy (EME) methods, as has been 
shown in [102] for finite-strain thermo-elasticity and in [106] for finite-strain thermo- 
viscoelasticity.’ 


> Of course, the GENERIC framework is not a prerequisite for the development of structure-preserving numerical 
methods for non-isothermal solid mechanics. In the context of coupled thermomechanical problems, alternative 
procedures have been proposed, such as in [72], [75], [76] and [73]. 


1. Introduction 


1.2. Objectives 


In this work, we aim to formulate, analyze, implement, and verify second-order structure- 
preserving numerical methods by employing the GENERIC structure. The present work 
can be divided into two parts: 


e The first part (Chapter 3) deals with structure-preserving numerical methods for 
molecular dynamical systems. It extends the scope of application for Energy- 
Momentum integrators to problems with periodic domains by making use of the 
discrete gradient operator [46], targeting the specific geometric structure of molec- 
ular dynamical systems. 


* The second part (Chapters 4 to 6) of this thesis adresses structure-preserving nu- 
merical methods for thermo-viscoelastic continuum dynamics. The focus hereby is 
on smooth irreversible effects, and therefore effects like plastic or damage transfor- 
mations are not considered. 


1.2.1. Energy-momentum conserving integration schemes for 
molecular dynamics 


Given the favorable properties of EM methods, it is remarkable that they have not received 
more attention in the field of molecular dynamics, at least when investigating the use 
of implicit time integration schemes for obtaining numerical solutions. The governing 
equations of the latter are essentially Hamiltonian and fit seamlessly in the framework 
developed since the 1970's for integrating these types of problems, while preserving 
the energy and the momenta. It would seem natural that integration schemes designed 
to preserve the main invariants of the motion would give accurate predictions of the 
thermodynamic averages, which are of interest in many practical and theoretical situations 
[107-111], but have rarely been studied [112]. Instead, molecular dynamics codes seem to 
favor the use of the explicit Verlet method or symplectic methods, mainly due to their 
smaller computational cost as compared with implicit schemes. Although popular explicit 
methods have desirable properties in terms of computational cost, accuracy, and geometry 
preservation, they lack energy conservation, a key invariant that is most important in 
the simulation of microcanonical ensembles. A thorough investigation of the accuracy 
of energy and momentum conserving schemes in capturing the statistical behavior of 
atomistic systems for long periods of simulation is lacking. Preliminary results [112] are 
promising, but much testing and validation is still required. 


The development and implementation of energy and momentum conserving algorithms in 
the context of molecular dynamics has three specific issues that affect the discretization of 
the equations and their analysis—issues that do not appear in their application to nonlinear 
solids, shells, rods, etc., or any of the other systems for which the use of these methods is 
widespread. The first critical issue is the treatment of periodic boundary conditions. These 
are almost invariably required for the study of average properties in particle systems [110, 
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111, 113], and demand a careful analysis, especially to ascertain whether they spoil the 
conserving properties of the method or not. Taking this into account requires to consider 
the geometry and topology of the periodic configuration space, and might also affect the 
accuracy of the integration scheme. 


The second issue that needs to be carefully dealt with is the use of conserving schemes in 
the context of three-body potentials. These functions are employed in modeling angle 
interactions in atomic bonding [114, 115], and their impact on the global behavior of 
some systems is so critical that it must to be accounted for. In fact, atomic systems with 
potentials of this type allow for large relative motions among the particles, and this is 
precisely the area where conserving schemes have shown their superiority with regard to 
other implicit integrators. 


The third aspect that requires a detailed analysis is the application of conserving schemes to 
mechanical systems in which the potential is based on cluster functionals [116-120]. These 
effective potentials are often required for the correct modeling of complex binding among 
metallic atoms and again require a careful study when used in combination with conserving 
schemes. While pair potentials of the Lennard-Jones type [121] have been employed 
together with EM conserving schemes [41], their formulation for cluster potentials needs 


to be specifically addressed. 


Here, we formulate energy and momentum conserving schemes for simulating the dynam- 
ics of atomic systems. Some of the methods discussed have already been employed in the 
literature, and we identify new ones. In all cases, we explain how the three critical issues 
identified before (i.e., periodic boundary conditions, three-body potentials and functional 
potentials) affect their formulation. To the author's knownledge, none of these have been 
previously studied. 


1.2.2. Structure-preserving integration of large-strain 
thermo-viscoelasticity in the framework of GENERIC 


The aforementioned GENERIC-based numerical methods are subject to two serious limita- 
tions: the first limitation is related to the use of the entropy density as the thermodynamic 
state variable, and the second limitation is due to the fact that the GENERIC formulation 
focuses on closed systems and thus does not account for boundary conditions. Of course, 
both points limit the applicability of numerical methods for the solution of initial boundary 
value problems in thermo-mechanics. 


The first limitation has been addressed in [122, 123] by using the absolute temperature as 
the thermodynamic state variable in the underlying GENERIC formulation. In the context 
of discrete systems, the free choice of the thermodynamic state variable has been recently 
addressed in [124]. The second limitation has been circumvented in [106] by applying the 
Lagrange multiplier method to enforce temperature boundary conditions. 
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We completely resolve the aforementioned limitations of GENERIC-based numerical meth- 
ods for finite-strain thermo-elasticity by proposing a novel GENERIC-based variational 
formulation that makes the free choice of the thermodynamic state variable possible. In 
particular, one may choose (i) the internal energy density, (ii) the entropy density, or 
(iii) the absolute temperature as the thermodynamic state variable. Moreover, boundary 
conditions are taken into account by applying a generalized GENERIC formulation for 
open systems. 


This novel, GENERIC-based variational formulation is discretized in time using the well- 
known mid-point rule. Depending on the choice of the thermodynamic state variable, 
partially structure-preserving schemes are obtained. For example, choosing the internal 
energy density as the thermodynamic state variable leads to an EM scheme. On the 
other hand, choosing the entropy density as the thermodynamic state variable yields a 
momentum-entropy (ME) scheme. However, despite their partially structure-preserving 
properties, all of the mid-point type schemes turn out to be prone to numerical instabilities. 
These observations led to the conjecture that only fully structure-preserving schemes 
guarantee numerical stability for dissipative systems in the same way that EM schemes 
do for Hamiltonian systems. 


In this context, the GENERIC-based weak form provides an ideal starting point for the 
development of EME schemes. First, the GENERIC-based weak form is discretized in space, 
resulting in a GENERIC-consistent space discretization. Then, the semi-discrete system is 
discretized in time by applying the partitioned discrete gradient operator in the sense 
of [125] leading to three different EME schemes. These EME schemes satisfy a specific 
Lyapunov-type stability estimate and thus do not exhibit any numerical instabilities. 
Hence, we develop novel EME schemes for large-strain thermoelasticity. 


Finally, we extend the applicability ofthis method to the realm of thermo-viscoelasticity. In 
particular, we aim at a material formulation of isotropic large-strain thermo-viscoelasticity 
which makes the free choice of the thermodynamic state variable possible. To this end, we 
build on previous work by [99], who laid the theoretical basis for the GENERIC description 
of thermo-viscoelasticity. 


1.3. Overview 


The thesis is organized into seven chapters. 
Each chapter can be read independently—therefore, repetitions are inevitable. 


Chapter 1 covers the motivation, the contemporary state of structure-preserving meth- 
ods, and the objectives of this work. 


6 A consistent notation throughout the present work would result in a tiresome reading experience, but 
notational shifts between the chapters are kept to a minimum depending on the focus therein. 
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Chapter2 introduces the GENERIC framework briefly. In Section 2.1 the framework 
is presented for finite dimensional systems. This includes its relation to Hamilton's 
equations of motion. Furthermore, we introduce the notation of GENERIC-consistent space 
discretization. Then, in Section 2.2, the structure of the building blocks of the GENERIC 
formalism is discussed when transformed to a different set of variables. 


Chapter 3 addresses the formulation and analysis of EM conserving time integration 
schemes in the context of particle dynamics, and in particular atomic systems. In Sec- 
tion 3.1, we review the basic topology of periodic systems in order to clearly define the 
distance function. Section 3.2 introduces particle dynamics, with special attention to its 
Hamiltonian structure. Next, in Section 3.3, conserving time integration schemes are 
presented for particle systems in periodic domains, restricted to those with pair potentials. 
These are extended to systems with angle potentials in Section 3.4 and to atomic systems 
with functional potentials in Section 3.5. 


Chapter 4 proposes a new GENERIC-based variational formulation for finite-strain 
thermo-elastodynamics. First, the GENERIC framework is dealt within bracket form for 
closed systems. Therefore, in Section 4.1, we start with finite-strain elastodynamics and 
subsequently present in Section 4.2 its extension to thermo-elastodynamics. Then, in 
Section 4.3, the transition to open systems is performed. In the following Section 4.4, the 
resulting GENERIC-based weak form of the initial boundary value problem (IBVP) at hand 
is then discretized in time and space. Representative numerical examples are presented in 
Section 4.5. 


Chapter5 designs novel EME methods for finite-strain thermo-elastodynamics. Section 
5.1 briefly recapitulates the variational formulation of large-strain thermo-elasticity. Then, 
in Section 5.2, the GENERIC-based weak form is discretized in space, resulting in a 
GENERIC-consistent space discretization. In addition to that, the main balance laws to be 
preserved under discretization are outlined. In Section 5.3, the semi-discrete system is 
further discretized in time, leading to three alternative EME schemes. Section 5.4 contains 
representative numerical examples that confirm both the structure-preserving features and 
the enhanced numerical stability of the newly developed EME schemes when compared 
to the standard time integration schemes developed in Chapter 4. 


Chapter6 introduces in Section 6.1 the bracket form of the GENERIC formalism for 
thermo-viscoelasticity with heat conduction. In particular, in Section 6.2, a new material 
version of the inelastic dissipative bracket is proposed which allows for the free choice 
of the thermodynamic state variable. In Section 6.3, a mixed finite element approach is 
proposed that yields a GENERIC-consistent semi-discrete form of the evolution equations. 


1. Introduction 


The mid-point type discretization in time is dealt with in Section 6.4, leading to alterna- 
tive partially structure-preserving schemes. Numerical investigations are presented in 
Section 6.5. 


Chapter 7 summarizes the findings. The thesis closes by drawing conclusions and 
giving an outlook on interesting research topics that follow the direction of this work. 
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2. A brief introduction to the 
GENERIC framework 


The general equation for non-equilibrium reversible-irreversible coupling framework 
(GENERIC) is a double-generator formalism that expresses thermomechanical models 
of dissipative materials (or generalized standard materials) in a unifying framework. Its 
formulation is based on an additive decomposition of the time-evolution equations into a 
reversible part and a dissipative part. While the reversible part is generated by the total 
energy of the system, the irreversible part is generated by the total entropy. 


The GENERIC framework is typically presented in two alternative forms: an operator 
version (i) in which the entries of the operator matrices are either generalized functions 
or differential operators, and (ii) a bracket version. Moreover, the GENERIC formulation 
typically focuses on closed systems. That is, neither thermal nor mechanical interactions 
with the surrounding ofthe system are considered. This Chapter summarizes the GENERIC 
framework for finite dimensional systems, which will be repeatedly referred to throughout 
the following chapters. For a comprehensive account of the GENERIC framework see, e.g., 
[94]. 


The fundamental concept of the GENERIC framework for finite dimensional systems also 
applies (with corresponding modifications) to infinitesimal dimensional systems, which 
will be considered in Chapters 4 to 6. In this connection the notion "GENERIC consistent 
space discretization" is introduced, which facilitates the design of energy-momentum- 
entropy (EME) consistent schemes. For now we will consider the case of isolated systems, 
for which the GENERIC framework was originally developed for. In the case of isolated 
(or closed) systems, boundary effects are considered to be irrelevant. 


2.1. GENERIC framework for finite dimensional 
systems 


The GENERIC framework hinges on an additive decomposition of the evolution equations 
into reversible and irreversible parts. Assuming that the state of a finite-dimensional 
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system is characterized by the state vector z(t) € IR", the time-evolution of the discrete 
system is described by 
dz 

qu Là,€ + MO,S , (2.1) 
which corresponds to the operator representation (i). In eq. (2.1) reversible processes 
are generated by the total energy €(z), while dissipative processes are generated by 
the total entropy S(z). In this connection, L(z) is the Poisson matrix, which must be 
skew-symmetric, whereas M(z) is the dissipative matrix, which must be symmetric and 
positive semi-definite. Two degeneration (or non-interaction) conditions have to hold. 
Namely, 


LO,S —0, (2.2) 
and 
Mo,é =0. (2.3) 


Using the GENERIC formulation (2.1) guarantees that the following two fundamental 
properties of a closed thermomechanical system are satisfied. Firstly, in accordance with 
the first law of thermodynamics, the total energy is conserved. That is, 

dE d 

— aE. E -0,€£ -LO,E + 0,£. MÓ,S = 0. (2.4) 
Note that on the right-hand side of the last equation, the skew-symmetry of the Poisson 
matrix has been used along with degeneration condition (2.3). Secondly, in compliance 
with the second law of thermodynamics, the total entropy should be a non-decreasing 
function of time. In fact, using (2.1) leads to 

dS dz 

— = 0,8: — =0,8 - LOZE + 0,8 -_Md,S > 0. (2.5) 
dt dt 
Here, degeneration condition (2.2) has been used, along with the positive semi-definiteness 
of the dissipative matrix. The properties of L and M can also be conveniently discussed 
in terms of two brackets 


{A, B} = 8,4: Lö,B, 


2.6 
[A, B] = 8,4. Mð,B . 2) 


Hereby are A(z) and B(z) arbitrary functions of the state variables. The brackets are 
termed Poisson and dissipative bracket, respectively. Using eq. (2.1) these brackets lead to 
the following evolution equation for an arbitrary function A 

dA 


ae {A,E}+[A,S], (2.7) 
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which corresponds to bracket version (ii) of the GENERIC framework. The antisymmetry 
property of the Poisson matrix can now be stated as 


{A,B} = —{B, A}. (2.8) 


Further, the Poisson bracket has to satisfy the product (or Leibniz) rule and the Jacobi 
identity. Those properties are well-known from the Poisson bracket of classical mechanics 
and express the essence of reversible dynamics, see e.g., [126]. 


The symmetry of the dissipative matrix M can be formulated in terms of the dissipative 
bracket as 


|A, B] = |B, A]. (2.9) 
The positive semi-definiteness can be expressed as 
[A, A] > 0. (2.10) 


Finally, the two degeneration conditions can be restated in terms of the bracket represen- 
tation 


(4,5) 20, 


=O. (2.11) 


2.1.1. Connection to Hamilton’s equations of motion 


If the system under consideration is purely mechanical and conservative, the choice 


z-qx xu uw coy p) 


where x“ is the a-th position, p^ the a-th momentum and N € Z4, results in a zero 
dissipative matrix M and a Poisson matrix L of the form 


L(z) = (2.12) 
—I 0 


where 0 is the zero and Lis the identity matrix with appropriate dimension depending on 
the system. By this means the state vector z has to be viewed as a column vector. The 
previous matrix is the well-known canonical symplectic matrix. In this case the evolution 
equation of the GENERIC coincide with Hamilton’s equation of motion, which will be 
considered in Chapter 3. 


13 


2. A brief introduction to the GENERIC framework 


2.1.2. GENERIC-consistent space discretization 


If an infinite dimensional system at hand is discretized in space in such a way that 
the GENERIC structure for finite dimensional systems is preserved, then we speak of a 
GENERIC-consistent space discretization. Such a method is devised in Chapters 5 and 6 of 
the present work. 


2.2. Transformation of variables 


By reasons of taste or convenience for the choice of the specific variables, it is important to 
address the relation between the building blocks of the GENERIC for different sets of state 
variables. We therefore consider a one-to-one transformation z +> z’. As for the scalar 
energy and entropy generators € and S we assume a simple transformation behavior 


£(z) = £(z(z)), 


2.13 
SE) = Sta). nd 
Applying a transformation of variables in eq. (2.1) using eq. (2.13) results in the following 
transformation laws for the Poisson and the dissipative matrix 


L'(z/) = 9,z - L(z(z’)) - 0,2! 


M'(z^) = 8 - M(z(z')) -ð d 
Considering that 
A'(z/) = A(z(z))) , 
Be) = Bis) a 
holds then the tanisformation laws, O44) lead to 
(AB) = Oy Al -L'OyB' = (sA: LOB) lou = {A B} C10) 
and 
LA^, BY] = Oy, A - M'ÓB' = (zA: MaB) aaa) = A Blare (217) 


In other terms: the Poisson bracket and the dissipative bracket are independent of the 
particular choice of state variables. 


The time evolution of the system is invariant with respect to the choice of a different 
set of variables. However, the building blocks of the GENERIC formalism, the matrices 
L, M and the gradients 0,€, 04S may have a simpler form, which can be of advantage 
for the formulation of a time-stepping integration scheme. Further, as an one-to-one 
transformation preserves the brackets, the transformations can be very useful when 
considering the verification of the Jacobi identity (see, e.g., [94]). 
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3. Energy-momentum conserving 
integration schemes for molecular 
dynamics! 


In this chapter we study the dynamic motion of systems of particles. This type of problems 
is of great importance for the simulation of matter at the atomic scale and a very large 
body of references study details pertaining to their numerical solution and the information 
that can be extracted from these simulations. 


Such systems can be considered as a prototype of a classical mechanical system: N 
particles in a periodic domain. The state vector of such a system assumes the form 


zc x. nu cx Bs. ED 


to be introduced in the subsequent. Hereby the state vector has to be viewed as a column 
vector. As only reversible processes are present, the GENERIC evolution eq. (2.1) boils 
down to Hamilton's equations of motion 

dz 

qu LE, (3.1) 
where L is the symplectic matrix eq. (2.12), or equivalently the equations of motion in 
Poisson bracket form 


dA 

—={A,E}. 3.2 

x Us (82) 
Therefore, solely the total energy E (which is the Hamiltonian 71 in this case) in connection 
with the Poisson bracket {-,-} describes the evolution of such systems, where the Poisson 


bracket is given by 


N 
{A,B} = V (kA: 054 B — Ops A xa B) , (3.3) 
a=1 
which is the well known Poisson bracket of classical mechanics. Depending on the specific 
structure of the total energy, different interatomic interactions can be modelled. Further, 


1 This chapter is based on [3] 


15 


3. Energy-momentum conserving integration schemes for molecular dynamics 


the GENERIC/Hamiltonian framework guides the construction of structure-preserving 
numerical methods, as we will show in the following sections. 


We address the formulation and analysis of energy and momentum conserving time 
integration schemes and identify three critical aspects of these models that demand a 
careful analysis when discretized: (i) the treatment of periodic boundary conditions, (ii) the 
formulation of approximations of systems with three-body interaction forces, and (iii) their 
extension to atomic systems with functional potentials. These issues, and in particular 
their interplay with energy-momentum (EM) integrators, are studied in detail. Novel 
expressions for these time integration schemes are proposed and numerical examples are 
given to illustrate their performance. 


3.1. The topology of periodic domains 


Many particle systems of interest are formulated in periodic domains. These allow to study 
large systems by only discretizing representative volumes, much smaller in size, while 
hopefully not losing too much information. In this section we gather some topological 
and geometrical facts of periodic domains that will be necessary to analyze numerical 
methods. 


We start by considering a periodic three-dimensional box B of side length L, noting that 
all the results are applicable to systems in one and two dimensions, with the corresponding 
modifications. This box is isomorphic to the torus T? (see e.g., [108, 110]), which itself 
can be identified with the product manifold S! x $1 x $+. Hence, each point € € T? can 
be uniquely characterised by three angles (a, 6, y) and the complete manifold is covered 
by a single chart. 


Numerical methods defined on T? pose difficulties that can be alleviated by mapping this 
set into a more convenient one. For that, let us first define the following equivalence 
relation on R°: two points x, y € IR? are defined to be equivalent, and indicated as x ~ y, 
if there exists a triplet of integers z € Z? such that x = y + Lz. Using this equivalence 
relation, we can define the quotient space P :— R?/Z that is homeomorphic to the torus. 
In what follows the equivalent class of a point x € IR? will be denoted as [x] € P. 


When dealing with systems of particles in periodic domains one has to choose one of 
the two homeomorphic descriptions described above, namely, the torus and P. From the 
computational point of view, employing the latter has many advantages. The first one is 
that given a distance on R°, this quotient space naturally inherits a distance, and thus a 
topology. In terms of the standard Euclidean distance d(-,-) : R? x R? > Rt U {0} we 
can define dr(-,-) : P x P++ R* U (0) by the relation 


dr([x],[y])= inf d(x,y). (3.4) 
x€ [x], yey] 
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Abusing slightly the notation, from this point on we will write d (x, y) instead of 
dr ([x], [y]). The second advantage of such a choice is that, for each equivalent class 
[x], there exists a unique point x € [x] n [- L/2, L/2)? that serves as identifier of the 
whole class which, in practical terms, implies that all operations need to be performed 
as with standard points in a cubic box. This identifier can be found using a projection 
operator 


7: ROB, (3.5) 


defined as 


; + L/2 
Tj = T(X)j = £j — floor (225) L, (3.6) 


where z;, j = 1,2,3 denote the Cartesian coordinates of the point x and floor: R — Z 
is the function that gives the largest integer smaller than or equal to a given real number. 
See Fig. 3.1 for an illustration of the projection operator. 


81,2 -L/2 0 L9 3L/2 


Figure 3.1.: Graph of the projection operator 7 restricted to one of the three coordinates of Euclidean space. 


The projection map 7r determines some of the properties and limitations of the numerical 
methods employed on systems with periodic boundary conditions, and it is helpful to 
summarize some of its properties: 


1. v is a nonlinear, surjective, projection, that is m o m = 7. 


2. The point 7 (x) is the closest one to the origin among all points in [x], that is, 


a(x) — arg inf d(x,0). (3.7) 


x€[x] 


3. In general, for arbitrary x, y € R3, 


dr(x.y) = |r (x — y)| # Im) — my) - (3.8) 
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4. The map 7 is C? except on the planes x; = L/2-- k;L, with k; € Zandi = 1, 2,3, 
where it is discontinuous. Away from these planes, the gradient O7 is the identity 
I: R? > R?. A one-dimensional illustration of the gradient Ö, is given in 
Fig. 32. 


—3L/2 —L/2 0 L/2 31/2 


aj 


Figure 3.2.: Graph of the gradient of the projection operator 0x7 restricted to one of the three coordinates of 
Euclidean space. 
5. For any x € R? 


a(x) = —m(-x). (3.9) 


6. For any x,y € R? and a € IR? 
dr (x + a,y F a) = dr (x, y) = (3.10) 
However, in general, if Q € SO(3), 


dr (x. y) # dr(Qx, Qy) . (3.11) 


Hence, in contrast with the Euclidean distance, the function dr (-, -) is not invariant 
under the action of the special orthogonal group. 


In this section we provide the basic ingredients that describe the dynamics of particulate 
systems. This dynamical system is governed by Hamilton's equations of motion, and most 
of its complexity comes from the particle interactions, as given by the potential energy. 
Here we present a fairly general class of potentials that will be examined more carefully 
in Sections 3.3, 3.4, and 3.5. 
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3.2. Particle dynamics: basic description 


3.2.1. System description 
Let us consider a system of N particles labeled a — 1,2,..., N moving inside a periodic 


box B. Let the mass of the a-th particle be denoted as m“, its position as x^, and its 
velocity as v^ = x^, where the dot indicates the derivative with respect to time. 


A system of particles such as the one introduced possesses a kinetic energy defined by 


1 
Tu sm lv"? : (3.12) 


Introducing the momentum p° = m“v“ of particle a the kinetic energy reads 


AA 
T= A ama P" , (3.13) 
and the potential energy is given by 
V -V(tL). (3.14) 


modeling the energetic interactions among all the particles. This potential energy is 
always a model of the true interatomic interaction and as such there exists a large number 
of simple effective potentials that have proven their value in different contexts (gases, 
fluids, organic molecules, metals, etc.). 


3.2.2. Equations of motion 


The dynamics of systems of particles may be expressed in the GENERIC format which is 
governed by 


— =L0,€, (3.15) 
where the total energy € = T + V is given as the sum of kinetic and potential energy. 
As no irreversible effects are present the GENERIC framework coincide with Hamilton’s 


equation of motion (see Section 2.1.1). Eq. (3.15) gives rise to the following evolution 
equation for the state variables 


(3.16) 
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The gradient of the potential energy is related to the forces acting on the particles, and 
we define 


f* = VL), (3.17) 


to be the resultant of all forces applied on particle a. 


These standard equations need to be carefully studied since the topology of B is not 
identical to that of Euclidean space and the notion of derivative has to be re-examined. For 
the moment being, let us assume that this object is well-defined, deferring until Section 3.3 
a more detailed inspection. 


3.2.3. Conserved quantities 


Eqs. (3.16) describe the motion of systems of particles that often possess first integrals, that 
is, conserved quantities along their trajectories. These quantities are of great relevance 
to understand the qualitative dynamics of the system, to develop controls, etc. These 
momentum maps are related to the symmetries of the equations, according to Noether's 
theorem (see, e.g., [126]). We review them very briefly, since they have a direct impact in 
the formulation of conserving schemes. 


We consider only potential energies with translational invariant, that is, functions V such 
that for every vector c € R? satisfy 


Vita) = V(x* +e} ala) - (3.18) 


Differentiating both sides of this equation with respect to c and setting later c — 0 we 
obtain the relation 


N 
o=) 3e VU) e. (3.19) 


This invariance condition is related to the conservation of the linear momentum of the 
system, defined as 


L- p°. (3.20) 


The time derivative of this quantity follows from the definition of particle momentum 
and eq. (3.19): 


N N 
b= pc. (3.21) 
a=1 a=1 


20 


3.2. Particle dynamics: basic description 


Systems of particles moving in the Euclidean space IR", with n — 2 or 3, often conserve 
angular momentum. This is a consequence of the rotational invariance of the potential 
energy. However, systems defined on a periodic domain do not preserve it, in general 
(see, e.g., [127, 128]). One way to explain this loss of symmetry is by noting that the 
projection (3.6) is not rotationally invariant (see eq. (3.11)) and thus when used in the 
definition of the potential energy, it spoils the invariance of the whole system. 


The total energy of the system is given as the sum of kinetic energy T' and potential 
energy V. The time derivative of the energy can be evaluated using the equations of 
motion (3.16), giving 


N N N 
£e mw NX oy Pv y fk 0, (3.22) 
a=1 a=1 a=1 


proving that the total energy must be a first integral of the motion. Given the relevance 
of the aforementioned conservation laws, numerical schemes have been proposed that 
attempt to preserve them. In particular, energy-momentum schemes, the ones under study 
in this Chapter, have been designed to integrate the equations of Hamiltonian systems 
while preserving both linear and angular momentum, in addition to the total energy. From 
the previous discussion, however, it follows that when dealing with periodic systems, one 
might focus on the preservation of linear momentum and energy, only. 


3.2.4. The hierarchical definition of the potential energy 


The potential energy of a system of particles is a function with the general form given in 
eq. (3.14) satisfying the invariance condition (3.18). A hierarchy of functions of growing 
complexity can be defined considering interactions involving an increasing number of 
particles. This is abstractly expressed as 


N N 
V-Vo >> Valea, xs) + $) Va(%a,%0, Xe) +--+, (3.23) 
a,b=1 a,b,c=1 
azb az bzc 


where V; is a function involving k—tuples of atoms and satisfying eq. (3.18). The for- 
mulation of accurate potentials is an active field of research and we limit our exposition 
to the most common types. The reader may consult standard references for a detailed 
motivation and derivation of other types (e.g., [120]). 


A convenient way of formulating potential functions that are translationally invariant is 
to include atomic interactions only via the distance between pairs of particles. In general, 
this would mean that the functions V; employed in eq. (3.23) must be of the form 


Vo(x4, xp) = Vo(d(xa, xb)) , 


- (3.24) 
V3(Xa, Xb, Xc) = V3(d(Xa, Xb), d(xo, Xc), d(Xc, Xa)) ; 
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and similarly for higher order terms. 


3.2.5. Dynamics in periodic domains 


Eqs. (3.16) define the motion of particles in periodic domains, but special care has to be 
taken with the definition of the potential energy and its derivative. 


With respect to the potential, we note that, when formulating the dynamics of particles in 
periodic domains, the distance function d(-, -) in eq. (3.24) should be replaced with the 
distance dT (:, -) defined in eq. (3.4). 


An aspect with important practical implications is that hierarchical potential functions of 
the form (3.23) are invariably defined employing a cut-off radius that effectively limits 
the number of particles that interact with those within that distance, in the sense of 
dz (^, -). Moreover, in order to avoid the singularities in the definition of the gradient 
of this distance, the cut-off radius is always chosen to be strictly smaller than L/2 (see 
Fig. 3.2). Equivalently, the dimension L of the periodic box must be selected larger than 
twice the cut-off radius. Under this condition, we observe that a collection of N particles 
in a periodic box B is a mathematical representation of an infinite domain consisting of 
boxes of dimension L x L x L that repeat themselves in the three directions of space. 


The formulation of energy and momentum conserving schemes in this kind of domains 
must take these two remarks into consideration, and we explore them in the following 
sections, starting from the simplest potential function possible. 


In this section we study the formulation of energy and momentum conserving algorithms 
for systems of particles in periodic domains where the potential energy includes only 
pairwise interactions. In terms of practical applications, only the simplest potentials 
belong to this class (for example, Lennard-Jones’). They are only accurate for modeling 
noble gases, but are very often employed for benchmarking and the study of numerical 
methods. 
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3.3. Energy-momentum methods for periodic systems 
with pairwise interactions 


3.3.1. Equations of motion 


We consider again a system of N particles in a periodic box B of side L with equations of 
motion (3.16) and an effective potential 


N 
1 Y a V5 
Ves » V (dr(x*,x°)) , (3.25) 
bZa 
where V : Rt U {0} > R. Using the definitions 


r^^ .— qm (x^ — x?) and pe = dp (x?, x^) :— |E®], (3.26) 


the forces (3.17) deriving from a pairwise potential can be written as 


N zab 
: ~ cab E 
Faser. with f”= VES (3.27) 
b=1 
azb 
For the following sections we further define 
r® =x- x’ and r® = d(x, x^) = |r®®|. (3.28) 


3.3.2. Time discretization 


We consider now the integration in time of the equations of motion (3.16) of a system in 
the periodic box B and effective potential (3.25). To approximate their solution we will 
employ implicit time stepping schemes that partition the integration interval [0, T] into 
disjoint subintervals [t,,,t,+41] with tn = n At, and At being the time step size, assumed 
to be constant to simplify the notation. In the algorithms defined below, we will use the 
notation x,, to denote the approximation to x(t,,), and similarly for the velocity. Moreover, 
the symbol fn+a will denote the convex combination (1 — o) fn + a fn+1 for any variable 
f and 0 € a € 1. 
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3.3.2.1. Mid-point scheme 


The canonical mid-point rule approximates the equations of motion (3.16) by the implicit 
formula 


Kr ni x a 
At i? Vn41/2 , 
a Vat d Vi y ab (3.29) 
m EN cm X, fip 
b=1 
a+b 


Here we introduced fo: the mid-point approximation of the force acting on particle a 
due to the presence of particle b, that is, 


zab 
b. Yrrfzab n+1/2 
fy .— V (P172) [ra^ | ’ (3.30) 
n+1/2 
with 
-ab b 
n41/2 z T (Xy 11/2 = X» 11/2) , T 
—ab =d a b ( E ) 
Ta+1/2 — T (X5 11/2: Xn41/2) : 
As in the continuous case, the condition 
b b 
fup = fin > (3.32) 


holds due to eq. (3.9). The properties of the mid-point rule are well known. For example, 
this method preserves the total liner momentum of the system, defined at an instant tn 
as 


Lig = Yo mvi. (3.33) 


L L M ve ve = 
nt1— un _ antl Yn _ ab _ 
mo xa P o ae oe N 999 
a=1 a,b=1 
a+b 


where we have employed eq. (3.32). 


3.3.2.2. Energy and momentum conserving discretization 


It is possible to construct a perturbation of the mid-point rule that, in addition to preserving 
the linear momentum of the system, preserves its total energy. The key ingredient of such 
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methods is the so called discrete gradient operator, an approximation to the gradient that 
guarantees the strict conservation of energy and momentum along the discrete trajectories 
generated by the integrator. 


Conserving integrators for problems in molecular dynamics have been studied since the 
1970's [40, 41, 46], although never for systems with periodic boundary conditions, with 
the exception of [112]. In all of these works, the conserving schemes are variations of the 
mid-point rule (3.29) of the form 


Xn+1 = x, a 
At may 
; k (3.35) 
a Vn+1 ^ Vn M 
mA = -Dye V, 


where Dxa is precisely the discrete gradient operator, an algorithmic approximation to 
the derivative 0,«. In analogy to expression (3.30), the discrete gradient defines a force 
contribution, to be specified later, such that 


N 
Do mU (3.36) 
b=1 
a+b 
If the following condition holds 
figo = fup + O(At?) , (3.37) 


the second order accuracy of the mid-point rule will be preserved. If we want the new 
method to preserve linear momentum we note from the previous section that it suffices 
that the pairwise forces f ds mimic the symmetry condition (3.32), that is, 
b b 
flgo = —fugo- (3.38) 
Any second order perturbation of £2 with this property will result in a second order 
accurate integrator that preserves linear momentum. The "classical" EM method is con- 
structed in such a way, and preserves, in addition to the total energy, the linear and 
angular momenta of the system, the latter being important in domains without periodic 
boundary conditions [46]. For problems in molecular dynamics interacting through pair 
potentials and posed on periodic domains, the “classical” EM method is based on the 
discrete gradient (3.36) with 
Vita — Vin? Yep HPR 
a4 (3.39) 


gab gab  , zab ^? 
Prati Tn Trl r 


ab _ gab ab „ab io 
f too mi Louis ;In41i) m 


where we have introduced the notation 


Ware, (3.40) 
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and where the appropriate limit must be taken in eq. (3.39) when |r2^ , — 74?| — 0. This 
form of the discrete gradient is frequently cited in the literature of integration algorithms 
(see, e.g., [112]) and is responsible for the conservation properties of the method, also in 
periodic domains. However, it is not the only possible form and, in fact, it can be shown 
that there are an infinite number of discrete gradients [129], some of which can be more 
easily extended to the periodic case. More precisely, a second class of EM schemes follows 


from a new definition of the algorithmic approximation of the pairwise forces: 


[rab [rab ab ab ab 
Va n YA x fup i (25 4 — rh ) 


ab ab „ab ._ gab 
flgo (r5 yt) = fur ! ab n, (3.41) 
re, — reb| 
n4-1 n 
where n is the normalized direction given by 
b ab 
T244—r 
n+1 n 

N= ab ab (3.42) 

Irz^a In | 


The expression (3.41) results from perturbing the force $5 in the direction that produces 
work and then imposing precisely that this work should coincide with the difference 
Ve = yo, an idea that has been previously employed in the literature in order to 
formulate EM methods [129]. The proposed definition corrects the pairwise force between 
the particles a and b with a “small” term in the direction n depending on unprojected 
relative positions. This is the result of the fact that the first equation in (3.35) is not posed 
in the quotient space but rather in the full R?. It might be argued that such a correction is 
nonphysical because it is not defined on the quotient space, where the problem is posed. 
While this is true, the velocity equation is not posed on this space from the outset, and 


the proposed correction results from this mismatch. 


The EM force given in eq. (3.41) is symmetric in a and b. Hence, the method defined by 
(3.35) and (3.41) preserves linear momentum. To show that the method indeed preserves 
exactly the total energy exactly, it suffices to take the dot product of the (3.35); with the 
left hand side of (3.35), and vice versa, and then add the result over all particles, that is, 


N N 
Y EN Kan. CY. Vn e Y me Vrai Ya 
At At "kie At 
a=1 a=1 
N N 3.43 
X MGE TEN ETE: LS 
At At algo At 
a=1 a,b=1 
azb 


Subtracting (3.43)2 from (3.43), gives 


Tan+ı — Tn ab: ne XR 

a ae (3.44) 
a,b=1 
azb 
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where the total discrete kinetic energy at time tn is given by 


N 
Toc Y. 
a=1 


In the spirit of eq. (3.41), further rewriting results in 


m^v$-:va. (3.45) 


Nl = 


Y po Xn41 7 Xn ime b 
[9 . n+ n p» f£% .(x9 —x?-— x? + x? 
2 algo At At E algo ( n+l n n+l n) 
a, b— a,b= 
azb a<b (3 46) 
A n 

ed 1 fe 2 ( ab an) 
ze At algo Tni ra). 

a,b=1 

a<b 


Therefore, a necessary and sufficient condition for energy conservation is that the follow- 
ing directionality condition is satisfied 


N N 
a a a 1 7 (za [^ (= 
5 fi : um = ra?) = 9 5 (Ver) e VG) x (3.47) 
a,b=1 a,b=1 
a<b a+b 


It is important to emphasize that the proof is based on the inner product between the 
algorithmic EM force and the unprojected relative position vectors. Finally, to show that 
the EM scheme (3.41) is indeed a second-order accurate method, it suffices to prove that 
the correction term in the definition (3.41) is of size O(At?). Making use of the relation 


ab T/a ria 
fyr = —O V (Ir 1 721) = pa V (ER 1721) ; (3.48) 
a Taylor series expansion around the point P j2 gives 
Vo Ve = fuera — rg) + O(AP). (3.49) 


Then, since the direction vector of the correction has size 


ry = ia 
—— —— =O(1 3.50 
Be] m > 


and [r?^ — r?^l is O(At), we conclude that the correction term is indeed O(At?). 


3.3.2.3. Time reversibility 


Time-reversible (or symmetric) integration schemes are often favored for the approxima- 
tion of Hamiltonian systems for two main reasons. First, the Hamiltonian flow itself is 


27 


3. Energy-momentum conserving integration schemes for molecular dynamics 


symmetric, so it is desirable that its numerical approximation also possesses this property. 
Second, symmetric numerical schemes are known to have several favorable properties [7], 
especially in long-term simulations. The class of EM integration schemes defined in this 
section have also this property. This is a direct consequence of the time reversibility of 
the algorithmic approximation of the pairwise forces, namely, 


b b b b b b 
Eie iEn (Eri) = en eee rn ) ’ (3.51) 


that is trivially satisfied by both (3.39) and (3.41). 


3.3.3. Interatomic potential 


For the following numerical examples we consider the well-known Lennard-Jones potential 


[121] with r = rè that is, 
m 12 6 
Hr) = (^ -(5)]. (352) 


where £ and c are material constants. 


3.3.3.1. Cut-off distance considerations 


In the Lennard-Jones potential, atomic interactions between distant particles are negligible. 
For this reason, a cut-off distance r. is often introduced beyond which the interaction is 
completely ignored (see, e.g., [108, 110]). However, simply trimming the Lennard-Jones 
potential beyond the cut-off distance leads to a discontinuity in this function at r = re that 
might affect the properties of the integration scheme. Since the derivative of the potential 
enters the equations of motion (3.16), this discontinuity precludes the computation of 
the interatomic force at r = r.. The discontinuity can be avoided, first, by shifting the 
potential function by the amount V (r..), leading to the shifted potential (SP) 


3.53 
0 ifr > re. 99) 


- V(r) — V(r fr<re., 
The derivative of this function at r = re is still not defined and neither is the force. To 
resolve this physical inconsistency one can introduce a shifted and linearly truncated 
potential (SF), which is equivalent to a shift in the force (see, e.g., [110, 128]), and given 


by 


(3.54) 


RN G -V (re) - (r= re) V’ (re) ifr < re, 
0 ifr > re. 
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One can further introduce a quadratic correction term that yields the shifted and quadrat- 
ically truncated potential (STF), which is equivalent to a shift and a linear truncation in 
the force, 


Vere ie k — V (re) - (r — re) V’ (re) — $(r = re) V" (re) ifr <re, (3.55) 


0 ifr > re. 


This potential is twice differentiable. Due to its higher smoothness, it is better suited for 
structure-preserving schemes than the standard potential since it eliminates numerical 
oscillations in the energy evolution that sometimes appear when employing non-smooth 
potentials. As mentioned in Section 3.2.5, it is important to stress out that the cut- 
off distance must not be greater than L/2 for consistency with the minimum image 
convention. Due to the quadratic term in the corrected potential, a cut-off radius of 
Te = 5a is suggested. 


3.3.4. Numerical evaluation 


All numerical examples are based on a set of dimensionless units. The implementation of 
periodic boundary conditions follows the details of Section 3.1 (see also [108, 110, 111, 
113]), and the order of magnitude of the time step sizes has been selected using standard 
criteria from molecular dynamics simulations (see, e.g., [130] and references therein.) 


3.3.4.1. Accuracy study 


In the first numerical example, we consider a two-dimensional box [-L/2, L/2]? with 
two particles. The initial positions and velocities are given, respectively, by 


x =(0,0)", vi =(0,0)", 

x? =(1.9,1)', v?=(5,0)', 
where the first particle is constrained to remain on the center of the box by holding fixed 
its degrees of freedom and only the second particle is allowed to move freely. For this 
simulation we used the Lennard-Jones potential with a simple spherical cut-off distance 


of re = 2.50. The numerical values of the remaining parameters of the simulation can be 
found in Table 3.1. 
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Material parameters € 20 Final configuration 
o 1 and trajectory 
Mass m“ 0.06 rn 
Side length L 12 
Newton tolerance - 10-6 "b 
Simulation duration T 0.8 
Reference time step size Atref 0.0005 
Time step size At 0.004,0.005,0.00625 
0.008,0.1,0.125,0.16 
0.02,0.025 


Table 3.1.: Accuracy study: Data used in the simulation 


Next we compare the “classical” EM (3.39), the newly defined EM (3.41), and the mid-point 
rule. For that, we compare the errors in the position and the linear momentum made by 
these three methods when using increasingly smaller time steps. In the two analyses, we 
use as reference the mid-point rule solution and define the errors, respectively, as 


a) lxt = xrl am _ IP“ — pr lle 


à = (3.56) 
° Ixzll á Ipzll2 


where 0$, with O € (x, p}, is the solution at time T calculated with the mid-point rule 
using the reference time step size At,-¢ and where [^ is the solution of the considered 
scheme at time T' for each time step size At. Fig. 3.3 confirm that all the schemes under 
consideration are second order accurate. 


1073 [ 


(MP) 


T 


1075 


—*— Midpoint 
—— EM 


—*— Midpoint 
—— EM 


—— EM classic —- EM classic 


Figure 3.3.: Accuracy study: Relative error in the position w.r.t. mid-point rule (left), Relative error in the linear 
momentum w.r.t. mid-point rule (right) 
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3.3.4.2. Energy consistency study 


The second numerical example investigates the energy conservation properties of the 
integrators described in Section 3.3. We consider 150 arbitrarily positioned particles inside 
a three-dimensional periodic box [— L/2, L/2]? such that the initial distance between the 
particles is greater than 2'/°o. Starting from rest, the kinetic energy of the system will 
rise until the system is in equilibrium due to the random order of the particles. For this 
simulation, we consider the shifted and quadratically truncated Lennard-Jones potential 
(STF) with a spherical cut-off distance of re = 50. Numerical values for the parameters of 
the simulation can be found in Table 3.2. 


Material parameters € 2 Initial configuration 
o 1 

Mass m“ 1 

Sidelength L 12 

Newton tolerance - 107° 

Simulation duration T 40 

Time steps At 0.08 


Table 3.2.: Energy consistency study: Data used in the simulation 


For this relatively large time step size, the mid-point rule introduces energy into the 
system leading to an energy blow-up and eventually to a termination of the simulation, 
indicated with a vertical line in Fig. 3.4 (left). 


19-10 
2 10 
1 
w 4 
0 
—1 
0 10 20 30 40 
t 
—— Midpoint —— EM 
—— EM —— EM classic 


— EM classic 


Figure 3.4.: Energy consistency study: Total energy (left), Total energy difference (right) 
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—— Midpoint 
—— EM 
— EM classic 


Figure 3.5.: Energy consistency study: Kinetic energy 


Both EM methods described, on the other hand, preserve the total energy up to machine 
precision for the whole duration of the simulation, see Fig. 3.4 (right). Hence, as long 
as the nonlinear iterative solver finds the step update, the solution of the EM will never 
blow up. Of course, if the time step size is fixed to value that is too large, the iterative 
solver will fail to solve the update equations and the method will not be able to proceed. 
Fig. 3.5 shows the evolution of the kinetic energy throughout the simulation. Since this 
energy is proportional to the temperature in the system, the figure reveals that only the 
EM methods can compute the evolution of the system until it reaches equilibrium, for the 
chosen time step size. 


We compare the newly proposed EM method with an explicit integrator, namely, the 
symplectic velocity-Verlet scheme provided by the molecular dynamics code LAMMPS? 
[131]. For this example we use the shifted and linearly truncated potential (SF) instead 
of the shifted and quadratically truncated potential (STF), since the former is available 
in the software package. Moreover, we extended the simulation duration to T' — 80. 
Keeping the time step size of the EM method Atgy = 0.08 constant, the time step size 
of the velocity-Verlet integrator At,y is reduced until the energy fluctuations become 
sufficiently small. As Fig. 3.6 reveals, the time step size employed for the explicit method 
is Aty = 0.0008, 100 times smaller than the time step size of the EM method. This is the 
time step size required to keep the relative energy fluctuations in the explicit solution 
approximately below 0.195. 


? https://lammps.sandia.gov 
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—— EM 
velocity-Verlet 
(LAMMPS) 


Figure 3.6.: Energy consistency study: Relative energy drift for Atem = 100At,y 


Regarding the computational cost per step, the proposed implicit integrator is naturally 
more expensive than the explicit scheme, but comparable to that of other second-order 
implicit methods like the mid-point rule. 


10° 


—— EM —— EM 


velocity-Verlet 
(LAMMPS) 


velocity-Verlet 
(LAMMPS) 


Figure 3.7.: Energy consistency study: Relative error in the position w.r.t. velocity-Verlet scheme (left), Relative 
error in the momentum w.r.t. velocity-Verlet scheme (right) 


In addition to the previous investigation, we perform an accuracy comparison between the 
EM and the velocity-Verlet methods. We use again the system described in Table 3.2 and 
obtain a reference solution with the explicit method and a small time step Atyer = 0.000001. 
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The displacement and momentum error measures of eq. (3.56) will be used to compare 
the accuracy of the two integrators for the time step sizes At = (0.0001, 0.0002, 0.00025, 
0.0004, 0.0005, 0.00625, 0.001, 0.0025, 0.004, 0.005, 0.01} and a final simulation time T = 4. 
As shown in Fig. 3.7, both schemes are second order accurate in the position and the linear 
momentum. 


3.4. Energy-momentum methods for periodic systems 
with three-body interactions 


In Section 3.3 we derived the equation of motion of a system of N particles moving inside 
a period box B where the interactions were based on pair potentials. For materials with 
strong covalent-bonding character, however, we further need to incorporate a bond-angle 
dependency to the effective potential. This can be archived by including three-body terms 
in the expression of the potential. 


3.4.1. Equation of motion 


We consider now a system of N particles in a periodic box B of side L with equations of 
motion (3.16) and an effective potential that depends only on the interactions among all 
triplets of particles. Such a potential must be of the form 


N 
1 P 
dar >, V (dr(x*, xh), dr (x^, x°), dr (x^, x*)) , (3.57) 


a,b,c=1 


abc 


where V : R+ U {0} x Rt U {0} x R* U {0} + R is a three-body potential between 
the a-th, b-th, and the c-th particle. The total kinetic energy of the system is given by 
eq. (3.13). The forces acting on the particles defined by eq. (3.17) can be written using the 
definitions in eq. (3.26) as 


N zab 
=)? with fp. (3.58) 


b=1 
a+b 
The strength of the force is now obtained by 


N 
b_ ) D(zab = zb 
yp = r Op V (r^ FT 2) i (3.59) 
c= 
azbzc 
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3.4.2. Time discretization 


Next, we consider the integration in time of the equations of motion (3.16) of a system in 
the periodic box B and effective potential (3.57) and employ the time integration strategy 
outlined in Section 3.3.2. 


3.4.2.1. Mid-point scheme 


The canonical mid-point rule approximates the equation of motion by the implicit formula 
(3.29), where the mid-point approximation of the force acting on the a-th particle in the 
direction of the b-th particle is given by 


red 
geb m" ab n41/2 
MP — PMP gao | | , 
n+1/2 
N (3.60) 
ab _ 5 V poe pac poe 
Yup = This n41/2:1n41/2: n41/2] > 


c=1 
azbZzc 


As in the continuous case, the weak law of action and reaction is satisfied and therefore the 
approximation preserves the total linear momentum of the system, see Section 3.3.2.1. 


3.4.2.2. Energy and momentum conserving discretization 


As in Sections 3.3.2.2, it is possible to construct a perturbation of the mid-point rule which, 
in addition to preserving the total linear momentum, preserves the total energy. Instead 
of the discrete gradient operator the partitioned discrete gradient operator [46] will now 
be employed, which is an approximation similar to the discrete gradient but is applicable 
to functions with more than one independent variable. To this end, let us rewrite the 
potential energy function (3.57) in terms of N(N — 1)/2 independent variables, e.g., 


V =V (F738, FAN p23,,,, pN- UN) = Ver), (3.61) 
where the double indexed set is given by 
(r^ = {7 | a,b (1,...,N), a < b}. (3.62) 


To simplify the definition of the partitioned discrete gradient, it proves useful to re-label 


the interatomic distances 7?^ using only their position in the array (77^). Therefore, a 


single indexed set is defined by 


(r^) = {7° | a € (L,..., N(N - 1/2) . (3.63) 
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Note that here an ordering of the (e pairs (a,b) has been established. For example, 
the map (a,b) — a could be chosen to be lexicographic (e.g., [132, pg. 43]). Then, the 
potential energy can be expressed, abusing slightly the notation, as 


V = ÖY. (3.64) 


For a potential function like this one, the discrete gradient operator is defined as 


N N 
a 1 a a 
Du. V=-Y fin =-), 5 (eh er = l (3.65) 
b=1 b=1 
a+b azb 


with the contributions 


ab _ gab N 4 5 
Enda = fme + m rab n, 
41 
a " , (3.66) 
a u za a a 
gab _ gab EN: 4 1) = Vale} fip (ar — Ya ) n 
n+l,n ^MP'' ab b , 
Irz^a mM 
for which we introduced the compact notation 
pra aa) _ (al -a—1l sa satl -N(N-1)/2 
Vieng l? ) = V (Thyes Tn T Tudor Tua ) ? (3 67) 
7a —QN _ ysl -a—1 <a zatl -N(N-1)/2 t 
ELA JH V (Tay nf Tn pts ln ys 


To show that the proposed integrator exactly preserves the total linear momentum, it 
suffices to follow the proof outlined in Section 3.3.2.2 and details are omitted. Similarly, to 
prove the energy conservation property, it is sufficient to show that 


N 1 N 
ab ab ab\ __ 3 (7 /zab -ac —bc [//5zab ac „bc 
) f ideo i (rapi — rh ) E 3! (UU s men mp) zu vr Ta; ln )) : 
a,b=1 ` a,b,c=1 
a<b azbZzc 


(3.68) 


A straightforward manipulation shows that the proposed method with algorithmic forces 
given by eq. (3.65) satisfies this condition. 


Remark 1. Many three-body potentials are expressed in terms of the bond angles 0*^^ at 
particle a, between the bonds ab and ac. Since the angle itself can be written in terms of 
the distances, that is, 


geac (3.69) 


2 pad pac poe = arcco (Fee)? + (ree)? a (Fr) 
er E 


the composition V og has again the structure of the potential (3.61) and thus the partitioned 
discrete gradient operator defined before can be employed without modifications. 
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Remark 2. In problems without periodic boundary conditions, using eq. (3.28), the EM 
method reads 


N N 
1 
Dev =-Ď f =- 25 (f$ + R l (3.70) 
er; iz 


foe H 
n,n+1 a _ pa ab ab ? 
M+ Th n+1 t fn (3 71) 
7a a 7a a ab ab à 
gab u v | ia 1) = Vi l £f. ) Prati T Ya 
n+1,n 7 a ab b: 
Thai — Tn Tni ct Tn 
where we used the compact notation 
a a 1 a-1l „a „at+l N(N-1)/2 
Vivaah ) -— V (Tas , Tn ST Tuna ^! ntl ) , (3 72) 
a a 1 a—l „a „a+1 N(N-1)/2 : 
n4-1 atr ) = AGER Thai Tn , Tn ) $ 


This EM method preserves the total angular momentum in addition to the total energy and 
the total linear momentum. It can be used, e.g., for the simulation of bonded three-body 
interactions between macromolecules. 


3.4.3. Interatomic potential 


For our numerical simulation we consider the Stillinger-Weber potential [115], which 
includes two- and three-body contributions 


V 1 p a ac c 

V Y ef” o) + 3i M eff" o T" [o i" [o) (3.73) 
a,b=1 a,b,c=1 
art agb 


The pair contribution is given by 


? D (Bí-*—$-P)g.(f) iff«a, im 


0 ifr>a, 


where the hats indicate the normalization of the distances by c. After composing with 
the law of cosines (3.69), the three-body contribution takes the form 


Í (P, pae. pec) a hP, $96. pec) $ hf, pee pr) ER h(f^*, pee par) : (3.75) 
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with 


pl)? 4 (p2)2_ (p38)? 2 
h(t, #2, #3) = (t DA ), (**) +4) ga(f1, 92) if?! c aandf? <a, 


pl 
0 otherwise . 
(3.76) 
Additionally, we introduce the functions 
^ ^ = 
f)-—exp(lf—a ; 
dut ) p ([ ] ) (8.77) 


gs (f, £?) = exp (af? -a|!-«[r-— a] !) . 


Here 0^«c = g(r, 72°, r?°) = arccos(—1/3) ~ 109.47? minimizes the function h given 
in eq. (3.76), that corresponds to the underlying diamond structure of silicon. From this 
reasoning it follows that the function h penalizes bond-angles which differ from the ones 
in this crystal structure. The parameters in the potential are A, D, a, £, o, q, p, A, and ». 


3.4.4. Numerical evaluation 


We consider now the numerical solution of systems of particles with the Stillinger-Weber 
effective potential. The pairwise contributions to the potential are discretized according 
to Section 3.3; the remaining three-body interactions are defined as in Section 3.4.2.2. 


3.4.4.1. Accuracy study 


We consider in this example a three-dimensional box [— 7/2, L/2]? filled with 5 particles. 
The initial configuration of the system has a particle at the center of the box and the 
rest of the particles form bonds with the first one with angle arccos(— 1/3). In addition, 
these four particles are at distance 1 from the center. Starting from rest, the motion of the 
system results from the non-vanishing interacting forces among particles away from the 
center. Further parameters of the simulation can be found in Table 3.3. 


38 


3.4. EM methods for period systems with three-body interactions 


Material parameters € 1 Initial configuration 
A 0 
B 0 
[o 1 
A 21 
p 0 
q 0 
a 1.8 
y 1.2 
Mass m“ 1 
Side length L 4 
Newton tolerance - 10-8 
Simulation duration T 0.7 
Reference 
time step size Atef 0.001 
Time step size At 0.01, 0.0125,0.016 


0.02,0.25,0.04,0.05,0.08 


Table 3.3.: Accuracy study: Data used in the simulation 


As in Section 3.3.4.1, we perform first an accuracy analysis of the EM integrator using the 
mid-point rule as a reference. This study is carried out using ten different time step sizes 
for both integrators and employing the error measures eq. (3.56). Fig. 3.8 reveals that both 
schemes are second order accurate. 


1078 107?  10-! 10° 
At 


—*— Midpoint —*— Midpoint 
—B— EM —H— EM 


Figure 3.8.: Accuracy study: Relative error in the position w.r.t. mid-point rule (left), Relative error in the linear 
momentum w.r.t. mid-point rule (right) 
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3.4.5. Energy consistency study 


We consider now 64 atoms inside a three-dimensional box [-L/2, L/2]?, initially arranged 
in a perfect diamond cubic lattice structure. 


Material parameters € 1 Initial configuration 
A 7.049556277 
B 0.6022245584 
o 1 
À 210 
p 4 
q 0 
a 2 
y 1.2 
Mass m^ 1 
Sidelength L 4 
Newton tolerance : 107°? 
Simulation duration T 8 
Time step size At 0.04 


Table 3.4.: Energy consistency study: Data used in the simulation 


This lattice will be disrupted during the simulation as we consider an initial velocity 
associated to each atom such that the total initial kinetic energy is approximately 768.22. 
For the chosen time step size, the mid-point rule clearly violates the conservation of the 
total energy, see Fig. 3.9, leading to an energy blow-up and finally to a termination of the 
simulation, indicated by the black line on the same Figure. 
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Figure 3.9.: Energy consistency study: Total energy 
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— Midpoint V, — EM V; 
Figure 3.10.: Energy consistency study: Potential energies using the MP rule (left), Potential energies using the 
EM method (right). V; refers to the i—body contribution to the potential 


The largest contribution to the algorithmic energy error is due to the mid-point approxi- 
mation of the forces generated from the three-body contribution of the Stillinger-Weber 
potential. It can be observed in Fig. 3.10 (left) that the potential energy of the two-body 
terms remains bounded, while the potential energy of the three-body terms increases 


unphysically causing the energy blow-up and, ultimately, the termination of the simula- 
tion. 


A 
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Figure 3.11.: Energy consistency study: Total energy difference (left), kinetic energy (right) 


In contrast, using the EM method both contributions to the potential energy of the 
system remain bounded, see Fig. 3.10 (right). As expected, the EM method preserves 
the total energy up to round-off errors, see Fig. 3.11 (left). Eventually, as illustrated in 
Fig. 3.11 (right), the evolution of the kinetic energy reveals that the EM solution reaches 
thermodynamic equilibrium for the chosen time step size, in contrast with the mid-point 
rule. 


3.5. Energy-momentum methods for periodic systems 
described by the embedded-atom method 


In addition to three-body potentials of the type described in Section 3.4, more elaborate 
and accurate potentials include multi-body effects through an environment-dependent 
variable, resulting in effective potential that are extremely common, for example, in the 
simulation of metals. One of the most popular interatomic potentials of this class is the 
one employed in the embedded-atom method (EAM) [116-118], whose use in the context 
of conserving schemes is analyzed in this section. The interatomic potential in this case is 
of the form 


1 x [7 a „b y M (=a 
V=5 2 Vreer FG. (3.78) 


The EAM potential consists of a pair potential contribution and electronic energies F (5°). 
The latter is due to the embedding of a-th atom in a homogeneous electron gas of density 
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p^ [120]. The background electron density function ° is a linear superposition of contri- 
butions from each neighbor atom such that the electronic energy F'(p^) of the a-th atom 
depends on the interatomic distance dr (x^, x^) to each neighbor. 


In the case of metals, the environment of each atom is a nearly uniform electron gas and 
therefore the embedded-atom approximation is reasonable. Since we already investigated 
pair potentials in Section 3.3, we now investigate the discretization of the forces due to 
the electronic energy, noting that the resultant forces in an EAM potential must include 
the forces due to the pair potential, as well. 


3.5.1. Equation of motion 


We consider a system of N particles in a periodic box B of side L with equations of 
motion (3.16) and an effective potential that depends only on the electronic energy of each 
particle, that is 


V=S OF |X o (dr(x^,x)|, (3.79) 


where V : Rt U {0} > R. Here, ga : R^ U {0} > R is a function of the relative 
interatomic distance which represents a spherical electron density field around the isolated 
a-th particle [120]. Using the definitions in eq. (3.26), the background energy density is 
then given by 


=X n (7). (3.80) 


The forces acting on the particles are defined by eq. (3.17) and have the standard form 


N zab 
f=) P5. with =f? = pL (3.81) 


b=1 
a+b 
where the strength of the force has now the structure 
p = E" (g*) gj (FP) + Fg, (F°) . (3.82) 


We observe that, for this potential contribution, the direction of the interatomic force 
depends only on the distance between the a-th and the b-th particle, while its strength is 
further determined by the background electron density at the a-th and b-th particle. 
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3.5.2. Time discretization 


We consider now the integration in time of the equations of motion (3.16) of a system in 
the periodic box B and effective potential (3.78) and employ the same time integration 
strategy as outlined in Section 3.3.2. 


3.5.2.1. Mid-point scheme 


The canonical mid-point rule approximates the equation of motion by the implicit formula 
(3.29), where the mid-point approximation of the force acting on the a-th particle in the 
direction of the b-th particle is given by 


nab 


b b ` n+1/2 
fip = MP Sab I , 
EVA (3.83) 


Php =F (Die) % avery zm E” (yp) 91 (soa) > 


and the mid-point approximation of the background energy density reads 


N 
z -ab 
Pue = y» (724 1/2) ` (3.84) 
b=1 
bza 
As in the continuous case, the weak law of action and reaction is satisfied and responsible 
for the conservation of total linear momentum of the system, see Section 3.3.2.1. 


3.5.2.2. Energy and momentum conserving discretization 


As illustrated in previous sections, it is possible to construct a perturbation of the mid- 
point rule which, in addition to preserving the total linear momentum, preserves the 
total energy. For that, we follow once more the steps presented in Section 3.3.2.2. The 
partitioned discrete gradient will be used again with a slightly different structure due 
to the nature of the embedded function. The partitioned discrete gradient assumes the 
form 


N N 
1 
b b b 
Dx. V = — > Page = — ` 3 (ae + ius ; (3.85) 
a,b=1 a,b=1 
azb a+b 


with the contributions 


Mab ab ab ab 
AFi n+ — fur ` (hpa T Ea) 


ab | gab n4 
Dod m FMP [pas c rab| B5 
n+l n (3 86) 
A Fob o foe x (r2^ EE rab) » 
fae "EA geb n+1,n MP n+1 n 
n+1,n = Imp + I5 


[roh s | 
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for which we introduced the compact notation 


ES ud ze POS cre sn = Far) 
EFL aita) le); (y 
AF an = Frain (Pai) — Fas (02) l 
+ FL Ga a 7 FUEL ER 2): 
and 
b— N 
nF) IH) + ga(77 P) - I, galas) » 
= =. (3.88) 
Prin?) = > gel) + gel) + I, golf“) 
= e=b+1 
The densities 55 ,, , , (r^^) and p?,, , ,, (r^^) are defined similarly. 


To show that the integrator exactly preserves the total momentum, it suffices to follow 
the proof in Section 3.3.2.2. The critical condition that a conserving scheme must satisfy 


reads now 


N N 

X fio ae = (FG) - F(gs)) (3.89) 
b=1 a=1 

a<b 


which is indeed satisfied by the proposed method. 


3.5.3. Interatomic potential 


We consider for our subsequent analysis the Lennard-Jones-Baskes (LJB) EAM model [118, 
133], which is the extension of the Lennard-Jones material model into the many-body 
regime of the EAM formalism [134]. The total potential energy of the LJB model is given 


by 


N N 

[/ 1 zab "EO 

ies > V(r )+ 5 F(p"). (3.90) 
a,b=1 a=1 
a+b 
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3. Energy-momentum conserving integration schemes for molecular dynamics 


The two body part has been introduced in Section 3.3.3. For the EAM contribution we 
further define 


y 
al 
DI 

Q 
= 

ll 


1 
564217" (In (p") - 1), 


1 N 
B= > Yoo (E), 
19-1 
bza 
I (Ir) zm D (-8 (c !|r^^| Er 1)) if [r^^ | X Te, 


0 otherwise , 


(3.91) 


where e, c, 3, A, re and Z; are material parameters. 


3.5.4. Numerical evaluation 


Since we already investigated pair potentials in Section 3.3, we focus on the EAM part of 
the LJB material model for the following numerical evaluations. 


3.5.4.1. Accuracy study 


A three-dimensional box [—L/2, L/2]? is now considered with 5 particles inside it. The 
particles form a unit body-centered cubic (BCC) cell of side 2, centered within the box. 
Starting from rest, the system builds up kinetic energy due to the fact that the particles in 
the exterior of the BCC crystal are not in equilibrium. Further parameters of the simulation 
can be found in Table 3.5. 
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Material parameters € 2 Initial configuration 
A 1 
o 1 e 
7 A N 
2 12 
Tc 3 e | 
Mass m^ 1 | ro | 
Side length L 6 | 22 v 
Newton tolerance - 10-6 u 
Simulation duration T 0.5 
Reference 
time step size Alger 0.0001 
Time step size At 0.0005,0.001,0.022 


0.0025,0.004,0.005 
0.01,0.02,0.025 


Table 3.5.: Accuracy study: Data used in the simulation 


Proceeding as in Section 3.3.4.1, we perform an accuracy analysis of the EM integrator 
using the mid-point rule as the reference solution. To study the convergence of the 
numerical solutions, we employ ten time steps of decreasing size and the error measures 
defined in (3.56). Fig. 3.12 confirms again that both the mid-point rule and the EM method 
are second order accurate approximations. 


1 1074 1078 107? 107! 
At 


—»— Midpoint —— Midpoint 
—B— EM —=— EM 


Figure 3.12.: Accuracy study: Relative error in the position w.r.t. mid-point rule (left), Relative error in the 
momentum w.r.t. mid-point rule (right) 
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3.5.5. Energy consistency study 


In this example we consider 108 atoms inside a three-dimensional box [—L/2, L/2]°, 
where the atoms are arranged in a perfect face-centered cubic (FCC) lattice structure. This 
FCC lattice is perturbed by the initial velocities of the atoms which are imparted in such a 
way that the total initial kinetic energy is approximately 3.251. 


Material parameters 


Mass 

Side length 

Newton tolerance 
Simulation duration 
Time step size 


At 


Initial configuration 


2 
1 
1 
4 


Table 3.6.: Energy consistency study: Data used in the simulation 


Further parameters of the simulation can be found in Table 3.6. 
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Figure 3.13.: Energy consistency study: Total energy (left), Total energy difference (right) 


For the time step selected, the mid-point rule clearly violates the conservation of total 
energy, see Fig. 3.13 (left). As in previous examples, this leads to an energy blow-up 
and finally to a termination of the simulation, indicated by the black line on the same 
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figure. In contrast, the EM method preserves the total energy up to round off errors, see 
Fig. 3.13 (right). Finally, and as in the previous examples, Fig. 3.14 shows that the EM, 
but not the mid-point rule, is able to evolve the system of particles up to thermodynamic 
equilibrium. 


10? 


+ 


— Midpoint 
—— EM 


Figure 3.14.: Energy consistency study: Kinetic energy 
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4. The GENERIC framework for 
large-strain thermo-elasticity! 


In this chapter we study the motion of large-strain thermo-elastic solids. Due the presence 
of heat conduction irreversible effects have to be considered. The continuum based 
description will resort to the GENERIC framework for infinitesmial dimensional systems, 
where the time-evolution of an arbitrary functional A can be expressed as 


d 
= ={A,E} + [A, S]. (4.1) 


This equation represents a 2-generator formalism in which the reversible part is generated 
by the total energy € of the system via the Poisson bracket {-, -}, while the irreversible part 
is generated by the total entropy S of the system via the dissipative bracket [-, -]. In contrast 
to Chapters 2 and 3 the brackets now imply an integration over continues labels rather 
then summation over discrete indices, and functional derivatives are considered rather 
then partial derivatives. A concise definition of these quantities will be presented in the 
sequel within the context of a Lagrangian description of large-strain thermo-elasticity. 


We propose a new variational formulation for finite-strain thermo-elastodynamics, which 
emanantes from the GENERIC formalism and makes the free choice of the thermodynamic 
state variable possible. In particular, one may choose the absolute temperature, the 
internal energy density or the entropy density as thermodynamic state variable. To solve 
initial boundary value problems, the GENERIC formalism is extended to open systems. 
The discretization in time makes use of the standard mid-point rule. Depending on the 
choice of the thermodynamic state variable, partially structure-preserving schemes are 
obtained. For example, choosing the internal energy as state variable yields a new energy- 
momentum consistent scheme. Thus the newly developed GENERIC-based weak form is 
particularly well suited for the design of (fully) structure-preserving methods as will be 
shown in Chapter 5. Furthermore, numerical investigations are presented which confirm 
the theoretical findings and shed light on the numerical stability of the newly developed 
schemes. 


The first part of the present chapter deals with the description of the GENERIC framework 
which was originally developed for isolated systems. In the case of isolated (or closed) 
systems boundary effects are considered to be irrelevant. 


1 This Chapter is based on [1] 
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4. The GENERIC framework for large-strain thermo-elasticity 


4.1. Thereversible part: Isothermal elastodynamics 


We start with the formulation of large-strain elastodynamics. Since elastic deformations 
are reversible the GENERIC evolution equation (4.1) reduces to the Hamiltonian form 
which is solely based on the Poisson bracket. 


p 


e2 


ej 


ea 


Figure 4.1.: Reference configuration B and deformed configuration (B, t) at time t. In Section 4.1 the focus is 
on the motion of isolated elastic solids. That is, external tractions acting on the boundary are disregarded. 


Consider a continuum body with material points X = Xe, in the reference config- 
uration B C R? (Fig. 4.1). Here and in the sequel the summation convention applies 
to repeated indices. Moreover, e4 denote the canonical base vectors in R3. Within the 
Lagrangian description of continuum mechanics the deformed configuration of the body 
at time t is characterized by the deformation map y : B x T +> R?, where Z = [0, T] is 
the time interval of interest. Accordingly, the placement at time t of the material point 
X € B is given by x = (X, t). Its velocity is defined by v = ¢, where a superposed dot 
denotes the material time derivative. 


The deformation gradient corresponds to the Jacobian of the deformation map 
F = ôx. (4.2) 


The Hamiltonian formulation of nonlinear elastodynamics is based on the introduction of 
the momentum conjugate to 2. To this end we consider the Lagrangian functional 


MORE f (p, F,@) dV , (43) 
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4.1. The reversible part: Isothermal elastodynamics 


which corresponds to the difference between the kinetic energy and the potential energy. 
The associated Lagrangian density function is given by 


, — 

Ilp, F, $) = spe: e - (W(F) bp), (4.4) 
where p : B ++ R4 is the mass density in the reference configuration, W : R?*? > R is 
a stored energy function which characterizes hyperelastic material behavior, and b : B +> 
IR? represents prescribed body forces, which are assumed to be dead loads. The conjugate 


momentum density is now defined by 


P = plp, F, p) = pp. (4.5) 


Moreover, the Hamiltonian density function results from a Legendre tranformation of the 
Lagrangian density function (4.4) with respect to the velocity leading to 


h(p,F,p)—-p:oe-l(e,F,o), (4.6) 


subject to relation (4.5). With regard to (4.4), the Hamiltonian density function reads 


las 
h(e,F.p)— 5p 'p:p* W(F)—-b:e. (4.7) 


The equations of motion can now be obtained by applying Hamilton's phase space principle. 
Accordingly, we consider the stationarity of the functional 


T 
Ap) = | [ ee - tto, E o) dVdt , (4.8) 
0 B 


with fixed endpoints in time (t = 0 and t = T > 0). The first stationary condition is given 
by 


Tq 
d 
DA(Go p) 3e — = | | (p: (6-59) - hp edp, F + eVöp.p)) avt 
B 


e=0 


0 
T 

zl (par - i [Mecho tiov dt 
0 NB B Ae 
0 


(4.9) 


53 


4. The GENERIC framework for large-strain thermo-elasticity 


where the variations dy : B +> IR? are subject to the endpoint conditions in time, dp, = 0 
and öpr = 0. The second integral on the right-hand side of (4.9) can be recast in the 
form 


d 
xx] Mer edp, +eVöp.p)av 
de B 


= i, (Oph - õp + Oph: Vow) dV 
e=0 B 


= / dp - (d.h — Div(Ogh)) dV 
B 

= 7 dp: ógchdV , 
B 


where 6,,h denotes the functional derivative (sometimes also called the Volterra derivative 
[135]) of the Hamiltonian density with respect to q. Note that in the above equation 
integration by parts has been applied along with the neglect of the boundary terms. As 
has been mentioned above, the neglect of boundary integrals is related to the fact that 
the focus in this section is on closed systems in which boundary effects are irrelevant. 
Performing integration by parts with respect to time on the first integral on the right-hand 
side of (4.9), the first stationary condition eventually yields 


Ey 
D ] e: set dVdt —0. (4.10) 
0 B 


The second stationary condition is given by 
d [T 
pay.p)-dp= 2 | f (p-cip) e — iG F, p cóp)) avar 
0 B 


" (4.11) 
x ] e (9-9 dV dt 

0 B 
=, 


Since the variations dp and dy are arbitrary (apart from the endpoint conditions on dp), 
the stationary conditions (4.10) and (4.11) give rise to the local form of the equations of 
motion 

p= óph , 


a (4.12) 
= —ôph , 


for all X € B and t € [0, T]. Taking into account the specific form of the Hamiltonian 
density function (4.7), the functional derivatives on the right-hand side of (4.12) are given 
by 

och = ph = Div(Oph) = —b- Div (OgW) ; 
and 


ôph pue p- pa 
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4.1. The reversible part: Isothermal elastodynamics 


Upon introduction of the first Piola-Kirchhoff stress tensor P = Of W, the local form of 
the equations of motion (4.12) can be recast as 


pp = Div (P) +b. 


This is the familiar Lagrangian form of the balance law for linear momentum (see, e.g., 
[136]). To translate the equations of motion (4.12) into the GENERIC bracket form (4.1), 
we now consider functionals of the form 


A(p,p) = ] 3e.» dv. (4.13) 


Differentiation with respect to time yields 
2d = nz F,p)dV 
= | (Gea: + dea: E +a: b) dV 
= [on ere) dV 
= f (Spa: dph — dpa: ph) dV =: { A, H} . 
Note that integration by parts has been applied. In this connection the boundary terms 


have been neglected again. The resulting equation 


dA 
— =A, H}, 4.14 
acd (4.14) 
holds for arbitrary functionals of the form (4.13) and is sometimes referred to as the 
equation of motion in Poisson bracket form (see, e.g., [24]). Note that this equation 
represents the reversible part of the GENERIC evolution equation (4.1). The Poisson 
bracket introduced above assumes the canonical form 


{A,B} =| (pa : dpb — dpa: dpb) dV , (4.15) 
B 


being the infinitesimal version of the classical discrete Poisson bracket eq. (3.3). Hereby 
are A and B two arbitrary functionals of the form (4.13). It can be easily observed that 
the canonical Poisson bracket (4.15) is skew-symmetric, that is 


{A,B} = - (B, A}. (4.16) 


Note that conservation of the Hamiltonian is an immediate consequence of property (4.16) 
implying (71, H} = 0, so that (4.14) yields dH. /dt = 0. Moreover, it is well-known that 
the Poisson bracket satisfies the Jacobi identity which is of paramount importance for the 
time-structure invariance of reversible dynamics (see, e.g., [94]). 
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4.2. Extension to dissipative systems: 
Thermo-elastodynamics 


Departing from the isothermal problem outlined above, we now aim at the inclusion of 
thermal effects. In particular, we focus on the dynamics of thermoelastic solids with heat 
conduction (see [136, Sec. 9.1] for a concise definition of a thermoelastic solid). We stress 
that we still focus on closed (or isolated) systems. 


In thermal solid mechanics the selection of the independent thermodynamic variable is a 
crucial point. From the viewpoint of constitutive modeling the temperature is a natural 
choice. The temperature field 0 : B x Z — R+ typically enters the definition of the 
Helmholtz free energy density 


V = ¥(F, 0). (4.17) 
Alternatively, the entropy density field 7 : B x Z +> R defined by 
n =T(F, 0) = —(F, 0) , (4.18) 


could be used as thermodynamic variable. Then, for example, the temperature can be 
expressed in terms of the entropy density by inverting the last equation leading to 


0 = O(F,n) . (4.19) 
Yet another viable alternative is to use the internal energy density field u : B x Z — R, 
which follows from a Legendre transformation of the free energy density with respect to 
the temperature. Accordingly, 


u = u(F,0) = 07(F,0) + v(F,0). (4.20) 


As usual, we assume that Ogu > 0 and ôs) > 0 for all (F,0) with det(F) > 0 and 
0 > 0. Using the internal energy density as independent thermodynamic variable, the 
temperature can be obtained by inverting the last equation leading to 0 = O(F, u). 
Making use of this relation, the entropy density (4.18) can be expressed in terms of the 


internal energy density via 

n — (E, O(F, u)) = f(F,u) . (421) 
Furthermore, inserting (4.19) into (4.20), the internal energy density can be expressed in 
terms of the entropy density, i.e. 


u=uf(F,n). (4.22) 


Asa further consequence of the Legendre transform (4.20) we get the following formula 
for the temperature in terms of the entropy density: 


6 = O(F, n) = d,u(F,n). (4.23) 


A clear exposition of further interrelationships between the alternative thermodynamic 
variables can be found in [100, Sec. 2.3]. 


56 


4.2. Extension to dissipative systems: Thermo-elastodynamics 


4.2.1. Free choice of the thermodynamic variable 


Guided by [100], in what follows, we allow for the free choice of the independent thermo- 
dynamic state variable by introducing the variable 7 : B x Z ++ R, which can be selected 
among the three aforementioned alternatives, i.e. 7 € (0,7), u}. In particular, in view 
of (4.19) and (4.22), the thermodynamic state variable 7 can be expressed in terms of the 
entropy density via 


T —T(E,n). (4.24) 


On the other hand, the entropy density can also be expressed in tems of 7, as can be 
observed from formulas (4.18) and (4.21). Accordingly, we introduce the map 7’ such 
that 


n= (F,7), (4.25) 


for r € {0,7, u}. The following relations will be needed in the sequel: 
(4.26) 


To proof these relations we consider the identity n = r/ (F, 7(F, )) which follows from 
(4.24) and (4.25). Since 0,1] = 1, we get 1 = 010,7, from which follows (4.26). Similarly, 
since Óp1) = 0, we have 0 = Opn’ + OL5/Og7, from which follows (4.26)2. Eventually, 
we recall a very important formula for the temperature that is given by (see also [94] or 
[100]) 


(4.27) 


Note that if the entropy density is chosen as thermodynamic variable (i.e. 7 = n), formula 
(4.23) is recovered. In the case 7 = 6, formula (4.27) yields the relationship 0 = 05u/0s7j 
which is an immediate consequence of the Legendre transform (4.20). 


4.2.2. Entropy-based brackets 


It is well-known that both the Poisson and the dissipative brackets are invariant under 
changes of the independent variables (see, e.g., [94] or Section 2.2, where this property 
has been verified for discrete systems). The Poisson bracket assumes a particularly simple 
form when the entropy density is chosen as the thermodynamic state variable. This is 
quite obvious since the reversible dynamics does not affect the entropy. Similarly, the 
dissipative bracket assumes its simplest form when the internal energy density is chosen 
as the thermodynamic state variable. This has been observed in [137] in the context of 
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a compressible, non-isothermal fluid. Moreover, this observation lies at the heart of the 
special form of GENERIC proposed in [100]. 

In the present work we start with the entropy-based formulation of thermo-elastodynamics 
and subsequently apply a transformation of variables to obtain the expressions for the 
brackets depending on the choice of the thermodynamic state variable 7 € {0,7, u}. In 
the entropy-based formulation we consider functionals of the form 


A= A(p,p,n) = ] ie. o mav. (4.28) 


Remarkably, in the entropy-based formulation the Poisson bracket retains the canonical 
form (4.15), i.e. 


{A,B} = [ (dpa: dpb — dpa - dpb) dV, (4.29) 


while the dissipative bracket featuring in the GENERIC evolution equation (4.1) is given 
by (see, e.g., [105]) 


[A, B] = / V (675,2) .O?KV (65,5) av. (4.30) 
B 


Here K = K(F, 7) is a positive semi-definite and symmetric second-order material 
conductivity tensor. Note that in the entropy-based formulation the temperature is given 
by (4.23). That is, in (4.30), O(F, y) = O,U(F, n). It can be easily observed that the 
dissipative bracket (4.30) is symmetric, i.e. 


[A, B] = [B, A], (4.31) 


and positive semi-definite, i.e. [A, A] > 0. 


4.2.3. Change of the thermodynamic variable 


We next aim at a uniform description of the Poisson bracket and the dissipative bracket 
in terms of the thermodynamic state variable r € {0,7,u}. To this end we consider 
functionals of the form 


DE I dermo (432) 
B 


and perform a change of the thermodynamic state variable. As has been outlined above 
the thermodynamic state variable 7 can be expressed in terms of the entropy density 
through relation (4.24), ie. 7 = T(F, n). Consequently, the two alternative expressions 
for the functional A given by (4.28) and (4.32) can be linked to obtain 


Alp, p.n) = ] ite Fo. m dV = ] iie Fo st av. 
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To get the functional derivative of A with respect to y we consider 


£ [aw +cdy,F + eVoy, p, T(E + eVdy,7)) dV T 
= J Ga «dad (Bea 4:0, dOpr) swap dy 
E f 5p O — Div(Opa! + a dV 
2 f TEE D 
B 


Here, again integration by parts and neglect of the boundary terms has been applied. 
Similarly, we obtain 


d 
& f ale, F, p, FE, + òn) av 
de B e=0 
= [ mosso gav 
5 
= f ônô,ădV . 
B 


To summarize, we have the following relationships 
dpa = dpa’ — Div(óza'Og7) , 
ópà = Opa 5 
bn @ = ôra' OnT , 


where, in view of (4.32), 6-a’ = 0,a’. The derivatives 0,,7 and ÖrT in the above equation 
can be substituted from (4.26) to obtain 


f 
bya = dpa’ + Div (oer) ; 


Orn 
Opa = gu. (4.33) 
~ Ó.a 
öna Duy 


These relationships can now be applied in the transformation of the brackets. Inserting 
(4.33), » into the entropy-based expression for the Poisson bracket (4.29), we arrive at the 
Poisson bracket formulated in terms of 7 € (0,7, u}: 


LA B = i (da - dpb’ — dpa’ - dyb’) dV 
B 


(4.34) 
. (ôa . (6,0! 
+f (Div (oer) - dpb! — dpa’ - Div (F260) dV . 
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T 7 (OM ul On 0,u! Orn! u 

n n O,u(F,g)  UF,n) 1 Onu(F, n) 0 Opu(F, n) 
9 —ÓgV(F, 0) 9 uF, 0) Ogri(F., 0) Ogu(F, 0) Ogr(F 0) Ogu(F, 0) 
u AEU) [AEN u  O4K(F.u) 1 Og fi (F, u) 0 


Table 4.1.: Specific relationships needed for the evaluation of the uniform expressions for the Poisson bracket 
(4.34) and the dissipative bracket (4.35), depending on the choice for the thermodynamic state variable 7 € 


{9,n, u}. 


Obviously, symmetry property (4.16) still holds. That is, (.A', B’} = —{B’, A’}. Substi- 
tuting (4.33)3 into expression (4.30) for the entropy-based dissipative bracket we get the 
dissipative bracket in terms of 7 € (0,9, u}: 


[A’, B’] = fy ($5) -(0'PKV (>) dV. (4.35) 


Note that in the last equation O' is given by formula (4.27) for the temperature and the 
material conductivity tensor is given by K = K(F, O/) = K'(F, 7). Moreover, symmetry 
property (4.31) is still satisfied being [.A', B'] = [B’, A’). 


Depending on the choice for 7 € (0,1, u}, specific expressions for the respective brackets 
can be obtained from (4.34) and (4.35). In this connection Table 4.1 provides a summary 
of specific formulas. With regard to Table 4.1 it can be easily concluded that the Poisson 
bracket (4.34) assumes its simplest form in the entropy-based formulation, while the 
dissipative bracket (4.35) assumes its simplest form in the formulation based on the 
internal energy density. 


Example: The Helmholtz free energy density (4.17) that will be used in the numerical 
investigations (see Section 4.5 and, for the following Chapters, Sections 5.4 and 6.5) is 
given by 


V(F,0) = vi (F) + vo(8) — (0 — 69)vs(J) , (4.36) 


where 


wi (F) = 5 (Fr -3- 2184-37 2°) + WrolJ) , 


Vs(J) = 38Wya(J) , 


Wal) = 233^ (og sy? + (T= 1) . 
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Here, J = det(F) is the determinant of the deformation gradient and ju, A are prescribed 
parameters, c > 0 is the specific heat capacity at constant deformation, is the coefficient 
of thermal expansion, and 0p is the reference temperature. We refer to [74, App. C] for a 
detailed investigation of the specific Helmholtz free energy density (4.36). It is worth noting 
that the Helmholtz free energy density (4.36) is consistent with the special case of linear 
thermo-elasticity (see Appendix A.1). It is now a straightforward exercise to calculate the 
quantities needed in Table 4.1. In particular, the temperature-based formulation (cf. 0-row 
in Table 4.1) yields 


Sipe (z.) IT 


T(E, 6) 
u(F, 0) = v4 (F)  c(0 — A) + 8ovs(J) , 
Ogr(F, 0) = P ; 
Ogu(F,0) — c, 
Op7(F, 0) = vS(J)cof(F) , 
F,0) = 


Opw1(F) + Boyz (J)cof(F) . 


Here, cof(F) = J(F)~T denotes the cofactor of the deformation gradient. The formulation 
based on the entropy density (cf. n-row in Table 4.1) gives 


n—wva(J) 


ü(F,n) = vi(F) + c (doe 


n—wv3a(J) 


— 2 + Oov3(J) , 


O(F, n) = oe” = ÖfF,n) 
Ogü(F, n) = gi (F) — (ÖF, n) — 6o) v4(J)eof(F) . 


Moreover, the formulation based on the internal energy density (cf. u-row in Table 4.1) 
leads to 


it dd wi) - AQ) + el) 


u + cho — V1 (F) - y = Da] Ji : 


[9 


AUF, u) = clog ( 


oi.) = ( 


AemE,u) = —[O(F,u)] pdr (E) + Qo C) cof) + v cof) 
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4.2.4. GENERIC evolution equation 


The GENERIC evolution equation (4.1) can now be considered in more detail by using 
the previously derived Poisson bracket (4.34) along with the dissipative bracket (4.35). 
Accordingly, we obtain 

dA’ 

tes {AE} + [A'S], (4.38) 
for functionals A’ of the form (4.32). Note that the evolution equation (4.38) represents 
a uniform description of thermo-elastodynamics in terms of the thermodynamic state 
variable r € {0, n, u}. In (4.38), the total energy of the system is given by the functional 


& (e. p,T) = I e (e, F, p, 7)dV , (4.39) 
B 
with associated density function 
1 
(pF, pT) = 5p '"p:ptu(F.r)-b.e. (4.40) 


Note that in the isothermal case the energy density (4.40) reduces to the Hamiltonian 
density (4.7). In addition to the total energy, the total entropy of the system acts as second 
generator in the GENERIC evolution equation (4.38) and is given by the functional 


S’(p,T) = Jene. (4.41) 


In the last equation, 7 : B x Z ++ Ris the entropy density which has been introduced in 
(4.25). For later use we mention that the following functional derivatives follow from the 
definition of the total energy (4.39) and the total entropy (4.41), respectively: 


dpe’ = Ope’ — DivOge' = —b — DivOgw' , 
tne = p = ptp, (4.42) 
6,e’ = 0,e' = 0,u' , 
and 
den = —Divdpn’ , 
dp’ =0, (4.43) 
às = Orn! . 
Now the left-hand side of (4.38) can be written as 
ZA = [teren dV 
= f (Apa! p+ Opa’: F + Opa! - p+ da't) dV (4.44) 


= ji (d,a'-p+dpa’-p+6,a’t) dV. 
B 
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4.2. Extension to dissipative systems: Thermo-elastodynamics 


Note that again integration by parts along with the neglect of boundary terms has been 
applied. The Poisson bracket on the right-hand side of (4.38) follows from (4.34) together 
with (4.40). Thus we obtain 


LAE / (5,a' - Ope’ — dpa’: dye’) dV 
B 


t 1 
+f (Div (oen) . Ape! — dpa - Div (Soen) dV. 


Similarly, the dissipative bracket on the right-hand side of (4.38) follows from inserting 
(4.41) into (4.35). Accordingly, 


[4,5] = i V (55) . (0)? KV (22) dV . (4.46) 


Next, we introduce the material heat flux vector Q : B x Z — R? through Q = Q'(F, 7), 
where 


(4.45) 


q = ev (227) 


oru’ 
= (O'?KV (a) en 
- -KVe'. 


Note that in the above equation use has been made ofrelationship (4.27) for the temperature 
field O'(F, 7). Now, the bracket (4.46) can be recast in the form 


ôa’ 

[A',S'] = | v( z) -Q'av. (4.48) 
B -u 

It is important to emphasize that the present brackets need to satisfy the fundamental 

degeneracy requirements of the GENERIC framework 


(A,S') 20, TE 
[AS E ees l 
The validity of these equations for arbitrary functionals A’ can be easily verified. The 
degeneracy requirements (4.49) ensure the thermodynamical consistency of the GENERIC 
evolution equation (4.38) in the sense that d£' /dt = 0 and dS’/dt = [S’, S'] > 0. Note 
that the last inequality follows from the above expression for the dissipative bracket 
together with the positive semi-definiteness of the material conductivity tensor K. Ac- 
cordingly, the total entropy is non-decreasing in the present case of isolated systems in 
which the net heating vanishes. 
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Taking into account relations (4.44), (4.45) and (4.48), the GENERIC evolution equation 
(4.38) can be recast in the form 


0 «f dpa’ - (p — Ope’) dV 
B 
1 . 1 : Oru 1 
+] dpa’: | D +ôpe + Div dern dV (4.50) 
B N 


y : da 6,a’ 
+f (di — Div (or) - Ope’ — V (5) : q) dV. 


This equation has to hold for arbitrary functionals A’ of the form (4.32). Consequently, 
the first two integrals in the last equation yield the local equations 


e-p!p, 


4.51 
p=b+DivP, er 


that must be satisfied for all X € B and t > 0. In (4.51)5, the first Piola-Kirchhoff stress 
tensor P — P'(F, 7) can be identified as 


P' = Ogu/ — Oder , (4.52) 


where O' is given by eq. (4.27). Concerning the third integral in (4.50), integration by 
parts and neglect of the boundary terms yields 


1 1 
fse (7 + — der : V(p7'p) + Diva’) dV =0. 
B on u 


à, 
Correspondingly, we arrive at the local form of the evolution of the thermodynamic 
variable 7 € (0,7, u}, given by 

Or Ou 
which has to hold for all X € B and t > 0. 


1 1 
qu Bo : V(p 1p) - ——DivQ’ , (4.53) 


Next, we summarize main ingredients of the present GENERIC-based formulation resulting 
from the specific choice of the thermodynamic state variable. Choosing 7 — 0 leads to 
the temperature-based formulation. With regard to (4.52), the first Piola-Kirchhoff stress 
tensor assumes the form 


P- Opu — 00g . 


Furthermore, (4.53) yields the temperature evolution equation given by 
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4.3. GENERIC for open systems 


The formulation in terms of the entropy density results from the choice 7 = n leading to 
the first Piola-Kirchhoff stress tensor 


P = pi. 

The evolution of the specific entropy is governed by 

; lpi à 

7 —-——DivQ. 

is) 

Choosing 7 = u yields the formulation based on the internal energy density with associated 
first Piola-Kirchhoff stress tensor 

P--00jj. 
Furthermore, the evolution of the internal energy density follows from 

ù = —O8gfj : V(p^! p) — DivQ. 

Remark 3. Apart from the dead loads b = —O,e’, non-potential material body forces per 
unit of reference volume, f : B x T +> R?, and a heat supply field per unit of reference 


volume, R : B x Z ++ R, can also be included into the present framework. To this end the 
supply term 


(A, F" = f (5 f an R) av, (4.54) 
B 


$ 
ru 


needs be appended to the right-hand side of the GENERIC evolution equation (4.38). 


4.3. GENERIC for open systems 


As pointed out in [94], there is no need to pay special attention to the boundary conditions 
if one is interested in the local field equations within an isolated system. Indeed, the 
majority of works dealing with the GENERIC formalism have been conducted in the 
context of closed systems. However, when it comes to numerical simulation, the boundary 
conditions are obviously of paramount importance. Consequently, we now aim at the 
inclusion of boundary effects (Fig. 4.2) into the GENERIC-based framework for finite-strain 
thermo-elastodynamics developed above. To this end, we build on ideas presented in 
[138] and apply them to the present framework for thermoelastic solids. Accordingly, the 
Poisson bracket (4.34) and the dissipative bracket (4.35) are now viewed as full brackets 
that can be decomposed according to 


(AB = [45 B ha + UB youn ? 


[4,5] = [A,B], + [A,B] ($33) 


bulk boun * 
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That is, the full brackets are split additively into bulk and boundary parts, respectively. 


€i 


ea 


Figure 4.2.: Reference configuration B with boundary OB and current configuration (B, t) at time t. External 
tractions t = PN act on the boundary of the current configuration. In addition to that, the heat flux across the 
current boundary is denoted by q = Q-N. 


This split can be accomplished by applying integration by parts to the full brackets. 
Following [138], the proper interpretation of the GENERIC evolution equation (4.38) is 
now given by 


dA’ 1 1 F 1 
= {A'E Jou + [A’, S"] 


a (4.56) 


bulk * 


The degeneracy requirements (4.49) for the full brackets transfer to the corresponding 
degeneracy requirements for the bulk brackets 


14 S Fea =0, 


ae" A (4.57) 


bulk — 
for arbitrary functionals A’. 


Invoking the symmetry properties of the full brackets, ie. {.A’, B’} = —{B’, A’} and 
[A’, B'] = [B’, A’), (4.55) leads to 


{A,B} puk = See puk u A B' roun T UB A ao) ? 


1 1 1 1 1 / / i (4.58) 
[AB jae = [BAT un (46: 8’) - [B',A'] 


boun boim) ` 
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4.3. GENERIC for open systems 


Substitution from (4.58) into (4.56) yields the GENERIC evolution equation for open 
systems 
dA’ 


Fag a —{E', dug iu [s', A’] 


dt bulk (4.59) 


= HA), E son A {EA M ng [4 d boun -— [S", A] "e 5 


In the sequel, the above equation will be used to formulate important balance laws and 
to derive the variational formulation of large-strain thermo-elastodynamics. It is worth 
mentioning that Remark 3 also applies to evolution equation (4.59). 


4.3.1. Specific brackets 


Next, we provide the specific brackets featuring in the GENERIC evolution equation (4.59) 
pertaining to thermoelastic solids. Proceeding along the lines of [138], we start from the 
full brackets derived above and apply integration by parts. In this connection, the goal 
is to get the derivatives (04a, dpa’, 0,.a') free of any spatial derivatives. To this end, we 
rewrite the full Poisson bracket (4.34) as 


/ 
t) f na tite (var Bi) 
B 


In 
, Ó.a/ 
p n (aed = t an) : pb dV . 


(4.60) 


Applying integration by parts to the second integral on the right-hand side of (4.60), and 
taking into account the additive decomposition (4.55)1, we obtain 


/ 
T (aoa! ee (ae! zs (av - oe 2220) av 


7] 

$ 

«f (oe — den) : VOSU AV, 
B T. 


(4.61) 


and 


$ 
(AB Lum ie pb- (dr - m 2 NdA. (4.62) 


Here, the vector N denotes the unit outward normal field on the boundary 06 of the 
reference configuration (Fig. 4.2). Similarly, applying integration by parts to the full 
dissipative bracket (4.35), yields the corresponding bulk bracket 


T 1 zb" 
[A,B] en z sa Div (te^ (Fa) dV , (4.63) 
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along with the boundary bracket 


0,a’ 
OB Ou! 


i] boun 


[A.B (0)? KV (5 =) dA. (4.64) 


oru’ 

It can be easily verified that the degeneracy conditions (4.57) are satisfied by the bulk 
brackets (4.61) and (4.63), respectively. Table 4.2 contains a summary of the specific 
brackets appearing in the GENERIC evolution equations (4.59) for open systems. The 
brackets in Table 4.2 can be easily obtained from expressions (4.61) through (4.64) by 
taking into account the generators in (4.39) and (4.41). 


Table 4.2.: Summary of the brackets featuring in the GENERIC evolution equations (4.59) for open systems. 


2a 


2 . 0,a’ 
EA a7 7 i, (» Opal pp: (aov en (ae Ec 
+ [y : VOpa’ dV , 
[S",.A'] sus =, Div o'yKV Ora dV 
x I ðu! ; 


E 0,a’ 
Ie eus L—— is p lp : * = Om > NdA , 
{EA oun E dpa’ P'NdA, 
Ts =| 5s UN. Q dA, 
o-a 
Is -f oN 0')*KV (35) dA 


4.3.2. Balance laws 


Next, we deduce the balance laws in material (or Lagrangian) form from the GENERIC 
evolution equation (4.59). We start with the total linear momentum of the continuum body 
defined by L — Jo p dV. In particular, we choose the density a’ = l, with le =€-p, 
where £ € R? is arbitrary and constant. This corresponds to the functional £^. = £ - L. 
Substituting Le for A’ in (4.59), the non-zero contributions on the right-hand side of 
(4.59) are given by TE pea = [ab : € dV and TER = - Jose : PNdA. 


Consequently, we obtain 


dL 
eee ( [avs [ praa) (4.65) 


68 


4.3. GENERIC for open systems 


Due to the arbitrariness of £ € IR?, (4.65) coincides with the balance law for linear 
momentum. Note that the parentheses on the right-hand side of (4.65) contain the resultant 
external loads applied to the continuum body (see also Fig. 4.2). 


The total angular momentum relative to the origin of the inertial frame is defined by 
J = fg p x pdV. We choose a’ = je, where jé = &: (p x p). Correspondingly, we 
substitute the functional Je = €-J for A’ in (4.59). The non-zero terms emanating 
from the right-hand side of (4.59) are given by (£', v = —J,b-(€ x ~) dV and 
TEN, Je NN = — fog(& x p): PN dA. Note that in the bulk bracket use has been made 


of the symmetry condition FPT = PFT. Consequently, the GENERIC evolution equation 
(4.59) yields 


dJ 
£- g E (fexpav+ [ox maa) ; (4.66) 


Since € € R? is arbitrary, the last equation corresponds to the balance of angular mo- 
mentum. Note that the parentheses on the right-hand side of (4.66) contain the resultant 
external torque about the origin (see also Fig. 4.2). 


Employing expression (4.39) for the total energy in the GENERIC evolution equation 
(4.59), we obtain 


dE’ 
Pr -AEE ha. + [SE] yan. 


23 (AE E Proun + {EE Proun d [£^] boun -— [SE] iux 
= TEE) 3r [S'E] u 1658 Nas E [£',5"] 
zZ ee T I 8 od) ? 


In the above equation use has been made of the decomposition of the full brackets in (4.55) 


boun 


along with the skew-symmetry of the full Poisson bracket and the degeneracy condition 
(4.49). Taking into account the specific brackets IE, p ME =— Jeg p p:PNdA 
and [s ST ss = 5 N: QdA, the above rate of change of the total energy can be 


recast in the form 


d 


se 'p-ptu av= [v-gav+ f (pp: PN- N-Q) dA. (467) 
dt B dB 


This is the balance of energy in material form which is in line with the first law of 


thermodynamics (see, for example, [139, Sec. 2.3]). 


Inserting expression (4.41) for the total entropy into the GENERIC evolution equation 
(4.59) yields 


ds’ 1 1 1 1 £ 1 1 f 
— = -{€',8'} + [SS] SE Y... IS] 


dt 
E [S',S'] E [S',S'] 


boun 


boun ? 
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where again use has been made of the decomposition of the full brackets in (4.55), degener- 
acy condition (4.49), and (S, Cus = 0. Since, as has been shown above, [S’, d > 0, 
and [S’, S’] boun = Jag (0) ! N - QdA the above equation implies 


dS’ 1 

— >- —N-QdA. 4.68 

dt ^ Jag 9’ Q ee) 
This corresponds to the Clausius-Duhem form of the second law of thermodynamics (see, 
for example, [136, Sec. 5]). 


4.3.3. Initial boundary value problem 


We next deal with the initial boundary value problem (IBVP) pertaining to large-strain 
thermo-elastodynamics. To this end, we decompose the boundary OB of the continuum 
body into a displacement boundary 0,,B, on which % = 2, and a traction boundary 0,B, 
on which PN = t, where @ and t are prescribed functions for t > 0 (Fig. 4.3). Moreover, 
pB U 0,8 = OB and pB N 0,8 = (). Similarly, for the thermal part we consider the 
subsets 0,B and 0,8, with the properties 0-6 U 0,8 = 0B and 0,61 0,B = 0 (Fig. 4.4). 
Here, the thermodynamic state variable is prescribed on 0, B, i.e. T = T, whereas the heat 


flux is prescribed on 0,B, i.e. Q N = @. 


The goal is now to determine the motion e : B x Z 5 R3, the linear momentum 
p: BxZ  R?, andthe thermodynamic state variabler : B x Z > R. The unknown fields 
are subject to initial conditions of the form y(-,0) = X, p(-,0) = pVo, and r(-,0) = rini 
in B. Here, V is a prescribed material velocity field and 7!"! is a prescribed field of the 
thermodynamic state variable r € (0,1, u}. The unknown fields are determined by the 
variational problem to be dealt with in the next section. 
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0,B 


€i 
e3 


Figure 4.3.: Mechanical part of the IBVP. Note that t — PIN denotes prescribed external Piola tractions acting 
on the current boundary expressed per unit area of the reference boundary ôs B. 


ea 


Figure 4.4.: Thermal part of the IBVP. Note that q = Q - N is the prescribed rate of heat transfer accross the 
current boundary expressed per unit area of the reference boundary 0B. 


4.3.3.1. Variational formulation 


To deduce the variational formulation of the present IBVP from the GENERIC evolution 
equation (4.59), we choose the specific density function a’ = w', where 


w = Wo: P +Wp: PWT. (4.69) 
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Here, wo, Wp : B > R? and w, : B > R are test functions that have to satisfy the 
boundary conditions w, = 0 and wp = 0 on 0,8, and w, = 0 on 0B. With the choice 
(4.69), the left-hand side of the GENERIC evolution equation (4.59) yields 


d 


aW = | nehme beret) dv. (4.70) 


On the other hand, the specific brackets on the right-hand side of (4.59) assume the form 
(cf. Table 4.2) 


(E, Wa =H f (» -wp + o'p: (we + Div (ow) —P': Yw») dv, 
B T 
[Sw] = -f + piv ((e^?kv ( 5.) ) av 
? bulk B (ey "NT ? 
{WE}, =O J aip (- ae) NdA, 
OB 71] 


-f wp: P'NdA, 
OB 


[W’,S’] boun = ye GN i Q' dA , 


1 = 
EAU -f. aN: (O)?KV (75) dA. 


Accordingly, the GENERIC evolution equation (4.59) gives rise to the equations 


o= fwe: (e - 7) dV, 


0= | (wp: $- b) +P! : Vwp) d - wp: P'NdA, 
B OB 


(4.71) 


along with 


; : Wr a ] ee Wr 
0 -f (w+ — Div (oer) -p p+ @ Div ((e^ev (22) dV 
—1 Wr / 1 1\2 Wr Wr / 
-| ——_O N-—N.(0’)*KV N. dA. 
a ot) ON (uL jen Ne) 


Integrating by parts twice, we arrive at an alternative representation of the last equation 
given by 


0 -f (w+ vtm : (Zn) V (5) . (Q)?KV (75) dV 


Wr 1 
+f gN Q'dA. 
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Taking into account expression (4.47) for the material heat flux, and the symmetry property 
K = KT, the above equation can be recast in the form 


= i E f MT d m jg 
o= | (vei vto 9): (Zn) v(zz) al oe 


Taking into account the above stated boundary conditions, (4.71) and (4.72) give rise to 
the following variational formulation of the present IBVP: 


(4.72) 


0= | we p-p'p)d, 


o=) (w (b — b) +P’: Vwp) dV — Wp-tdA, 
B 8,B 


= . EET NT / Wr Wau Wr _ 
o-f (vet vo p): (Eon) V (55) q) dV + 5 pau 


(4.73) 


These equations have to hold for all times t > 0 and for arbitrary test functions subject to 
the above mentioned boundary conditions. 


Remark 4. The variational formulation (4.73) can be easily extended to account for non- 
potential body forces and a heat supply field. To this end, the supply term introduced in 
Remark 3 just needs be included in the above derivation of the variational formulation. 


Remark 5. The majority of previously developed numerical methods for finite-strain 
thermomechanics rely on the use of the temperature as the thermodynamic state variable. 
Choosing T = 0 in the GENERIC-based weak form (4.73)3, leads to the variational equation 
for the temperature evolution given by 


Ò+ V(p^! (3 “oem v(#2)-a) avs f Me ATE 
Jv (p p) Og] u Og Q dB du 


The last equation can be recast in the form 


Woe En Oot -1 > woe — / Wo _ 
0 4 ] — —]- d ——qdA- 
Li (aa amp p) 2) v(3%) a} Vib Loom 0, 


or 


/ { vg (dad 4 0V (p^!) : dem) + Vus - Kvo} dV «f vegdA —0. (4.75) 
B ƏB 


In the last equation the modified test function vg = EXT has been introduced. In addition 
to that, temperature formula (4.27) has been used along with constitutive equation (4.47) 
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for the material heat flux vector. In essence, weak form (4.75) of the temperature evolution 
equation lies at the heart of numerous works on finite-strain thermomechanics such as, for 
example, [140-142]. Alternatively, differentiating (4.20) with respect to time, one arrives 
at the relation 


800 = 07 — 00g : F . 


Inserting the last equation into (4.75), and making use of the identity ~ = o !p, we 
obtain the entropy form of (4.75) given by 


[teen + Veo KVO} avs f veg dA =0. (4.76) 
B 0,B 


In essence, variational formulation (4.76) lies at the origin of alternative temperature-based 
numerical methods such as, for example, [143, 144]. It is worth mentioning that in (4.76), 
7] = 47(F, 0) is the time derivative of the “constitutive entropy" 7(F, 0). This is to be 
contrasted with the entropy-based formulation emanating from the GENERIC-based weak 
form (4.73)3. Specifically, choosing T = n, (4.73)3 yields 


f { wo -y (7,6) .6?KV (8) } dV +f , Fada =0, (477 
where © = Onu(F, 1). 


4.3.3.2. Balance laws revisited 


The balance laws dealt with in Section 4.3.2 can also be recovered from the variational 
formulation (4.73), by choosing specific test functions. For that purpose we confine our 
attention in this section to the pure Neumann problem (i.e. 0,B = 0,B = OB). 


Consider wp = ple, where Ue = €-p has been introduced in Section 4.3.2. Correspond- 
ingly, inserting wp = £ into (4.73)5, we directly recover the balance law (4.65) for linear 
momentum. 


Next, we choose wy = pjg and Wp = Opj¢, where je = E- (P x p) has been introduced 
in Section 4.3.2. Accordingly, inserting w, = p x € and wy = € x q into (4.73)1 5, 
and subsequently adding both equations, we recover the balance law (4.66) for angular 
momentum. To verify the balance of energy we choose for the test functions in (4.73) 


= ! Bret 
Wy = Ove’ , Wy, ——b, 
— 8. e "d 
Wp = p€, or Wp=p p. 
t 
w, =0,e , w, =0,u'. 
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Substituting these quantities into (4.73) and subsequently adding the three resulting 
equations, a straightforward calculation recovers the balance law (4.67) for energy. In this 
connection, the following relation is of importance: 


/ 
i (P : V(p ! p) + ôru't + V(p !p): (der )) dV 
B 


dr 
Oru’ ru’ 
m LE / T AE GE /l. —1 
Kin 
=P’ (4.78) 
-f (On! + Opu’ : F) dV 
B 


d / 
= — d x 
IE v 


Here, formula (4.52) for the first Piola-Kirchhoff stress tensor has been used, along with 
the identity ~ = o^! p. 


Concerning the balance of entropy, we insert w, = 0,1’ into (4.73)3, to obtain 


1: Si / 0,1 / Oct e 
0= 9,77 +V(p p):(Oggy)— V -Q’) aV + ——-Q:NdA 
B ru! 6) ru! 


B Oru 
u dr’ 1 s3 1 T 
-| (SE -v(g)-xv oO dV + ng RAM 
Here, use has been made of formula (4.27) for the temperature along with expression 


(4.47) for the material heat flux vector. Moreover, the identity ¢ = p^ ! p has again been 
taken into account. The above equation can be rewritten as 


ds’ 1 1 1 

"C _ | yl —).(e’P2KV dV '. NdA 4.79 

a [v (a) evsv(S)a- [oa Nat, — um 
————M————— 


20 


which recovers the second law of thermodynamics in the form (4.68). 


4.4. Discretization in time and space 


4.4.1. Discretization in time 


We first perform the discretization in time of the variational formulation (4.73). To this end 
we focus on a representative time interval [t„, tn+1] with corresponding time-step size 
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At = ta+ı — tn. The discrete approximations at times tn and t„+1 of the continuous vari- 
able (e), will be denoted by (e), and (¢),,+41, respectively. Moreover, the approximation 
of any state variable (e), at mid-point time ty = $(tn + tn41) is given by 


(n+i = 59s + (@)n41) ; 


Assume that the state variables o, p, : B > Rê andt, : BWR ™m € On, Mn; Un} 
are given. We aim at the determination of the corresponding state variables at time 
tn+1 (Yn41>Pn41>7™m+1) through the mid-point type discretization of the variational 
formulation (4.73) given by 


Wr 


Wr / > 
edP |- Q dV 4 Gee dA. 
i CETINA Duc t 9,8 Ont! |, n3 


+ 


(4.80) 


Note that depending on the choice for the thermodynamic state variable 7 € (0,1, u}, 
three alternative semi-discrete formulations result from (4.80). With regard to (4.52), 
the time-discrete version of the first Piola-Kirchhoff stress tensor in (4.80)2 assumes the 
form 


P| = PF.) 


ET 
= Lm — e 2 


A frame-indifferent thermoelastic formulation requires that the functions u’ and 7’ can be 
expressed as 


(4.81) 


ta+4 


(4.82) 
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where C — FTF is the right Cauchy-Green tensor written in terms of the deformation 
gradient F. Taking into account (4.82), we obtain 


Opu' |, Pe = 2F, | 10cv"(C], TN (T8341) ? 
US A ? (4.83) 
Orn l; eo 2F,10cr (Cl, 1 tet) ’ 
"+5 3 
together with 
OETAN Do ru” (C|, , Tag) 
eus HL : (4.84) 
IN lui = Orn (Chey (7241) , 


where C|, , denotes the right Cauchy-Green tensor evaluated in the mid-point configu- 
n+4 
ration. That is, 


Cl, , = Fapa Fari (4.85) 


n+4 
Moreover, in view of (4.47), the time-discrete version of the material heat flux vector 
/ B g E n 
Q o featuring in (4.80)3 is given by 
2 
£ 1 
Q l. =Q (Fati Tnt) 


(4.86) 


, 


where the time-discrete version of the material conductivity tensor K'| ps follows from 
n+4 
the relation 


K’ = K'(E siinga) - RC Tri) : (4.87) 


4.4.1.1. Semi-discrete balance laws 


We next investigate to what extent the balance laws outlined in Sections 4.3.2 and 4.3.3.2 
are inherited by the three alternative semi-discrete formulations developed in the last 
section. To this end we proceed along the lines of Section 4.3.3.2. We again confine our 
attention to the pure Neumann problem (i.e. 0, B = 0,B = OB). Inserting wp = € into 
(4.80)5, we directly obtain 


e cm cg (f oar f fn. 44) ’ 
At B oB ^ 
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where t, | P P' | be N. The last equation corresponds to the semi-discrete version of 
nili 


2 
the balance of linear momentum (4.65). Inserting Wọ = Pati X £and wp — £x ucl 
into (4.80); », and subsequently adding both equations, a straightforward calculation 


yields 


Judi — Ja pu 
e TS ce (fosa bae f esa x EnaA) 


which corresponds to the semi-discrete version of the balance of angular momentum 
(4.66). Accordingly, independent of the choice for 7 € (0,7, u}, the three alternative 
formulations at hand consistently reproduce the balance laws for linear and angular 
momentum for arbitrary time-steps At. 

. ma ME: =). =) . . . 
Inserting wp = —b, wp = p Pn+4; and w, = Owl, into the semi-discrete 


2 
formulation (4.80), and subsequently adding the three equations, we arrive at 


1 f —Ppa: nt — f 
je Pa+1 Phi Pu’ Po (o ds n+l n E Opu' |, m MS) dV 
B n+3 ni 


2 At At 
mu . Qui — Pn — ‘It Ber; er ) 
= f» Poe Pave f (p pupa uw 
where n+ 1= Q'| TEA N. The last equation can be viewed as semi-discrete version of 


the balance of energy (4.67) provided that the second integral on the left-hand side of the 
last equation is equal to 4 fau (Fn+1; 7541) — W (Fn, Ta) dV. However, for 7 € {0,n} 
this is generally not the case. In contrast to that, for 7 = u, the above equation reads 


1 —1 Pn+1 ` Pn+1 = Pn’ Pn N Un+1 — Un 
dV 4 dV 
5f At rn: 


Qui — Pn J —1 T = 
= | b. —— dV + i-t nr) dA, 
| At 8B (> Pe hun * 4) 


which indeed corresponds to the semi-discrete counterpart of the balance of energy 
(4.67). Accordingly, the formulation in terms of the internal energy density leads to an 
energy-momentum scheme. 


(4.88) 


Eventually, inserting w, = 0;7'|, , into (4.80)3, we obtain 
ni 


Tn — Tn Ea 
lo. co denle, :V(p LS) dV 


2 


S nl n2 / /i—1 nl m 
= jv (e E. 0 m K luu V (e E av- f © ne nti dA. 
————————— 


N: 
>0 
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The last equation does conform with the balance of entropy provided that the integral on 
the left-hand side of the last equation equals 


1 
ai] i ene) = (Fn, tea 


However, in general this is not true for 7 € (0, u}. In contrast to that, for 7 = n, the last 
equation yields 


Mn+1 — Mn 
—————— dV = 
av 


—1 


ml "T 
Vio dV — ji is Gn41 dA, 
EE heel 0B ep 7 


which corresponds to the semi-discrete counterpart of the second law of thermodynamics 
in the form (4.68). Accordingly, the formulation in terms of the entropy density leads to 
an momentum-entropy scheme. 


4.4.2. Discretization in space 


We apply standard isoparametric finite elements (see, for example, [145]) based on finite- 
dimensional approximations of the state variables at time t, given by 


PX) = M NYX) palt), PX) = M NX) palt), — (489) 
a=1 a=1 
and 
(X) = Y N*(X) Talt). (4.90) 
a=1 


Here, N^ : B > R denote the nodal shape functions and v, (t), p,(t) € R?, ta(t) ER 
are the respective nodal values at time t. Moreover, Nnoae denotes the total number 
of nodes in the finite element mesh. The standard (Bubnov) Galerkin approach relies 
on analogous approximations for the test functions in the variational equations (4.80), 
Wy, Wp and wz, denoted by wy, Wp” and wP. 


It is worth noting that the semi-discrete balance laws investigated in the last section still 
hold for the fully discrete finite element formulation as can be easily verified by proceeding 
along the lines of the last section. 
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4.5. Numerical investigations 


In this section, the three alternative mid-point type schemes newly developed in the 
present work are applied to representative numerical examples dealing with finite-strain 
thermo-elastodynamics. Depending on the specific choice of the thermodynamic state 
variable r € (0,9, u}, the following methods are applied: 


u | (EM), | Energy-Momentum consistent scheme 
n | (ME), | Momentum-Entropy consistent scheme 
0 | (M)o Momentum consistent scheme 


|| 
ll 


Concerning the constitutive equation (4.47) for the material heat flux vector, we assume 
thermally isotropic material, with material conductivity tensor given by 


K-kJC |. (4.91) 


Here, k is a prescribed coefficient of thermal conductivity and, as before, J = ,/det(C). 
Moreover, the Helmholtz free energy assumes the form (4.36). 


In the numerical investigations we shall focus on momentum maps associated with 
symmetries of the mechanical system at hand, and the balance laws associated with 
the two fundamental laws of thermodynamics. In this connection we also consider the 
functional 


L=E—- 0S. (4.92) 


According to [146] (cf. Appendix A.2), for certain types of environments, £ is a natural 
Lyapunov function and thus qualifies as estimate for the numerical stability of the schemes 
under consideration. Concerning the numerical solution of the resulting algebraic system 
of nonlinear equations we apply Newton's method together with the convergence criteria 
summarized in Appendix A.3. 
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4.5.1. Flying L-shaped block 


ta 
360.00 
2 
333.33 1 
05 246 810 
= 306.67 : 
i [s] 
e2 280.00 — Ff) 
€1 


Figure 4.5.: L-shaped block: Discretised block with initial temperature distribution and mechanical boundary 
conditions (left), load function over time (right) 


The first numerical example deals with the L-shaped block depicted in Fig. 4.5. The spatial 
discretization of the block relies on 117 tri-linear finite elements leading to 224 nodes. 
The initial temperature field is varying linearly over the height (x3 direction) of the block. 
In particular, at x3 = 0, the initial temperature is prescribed as 0,, while at x3 = h, the 
temperature is prescribed as 0). The whole block is assumed to be thermally insulated 
(q = 0 on 0,B = OD). Starting at rest, Piola traction vectors ta and t; are acting on two 
parts of the boundary surface of the block (Fig. 4.5). The external loads are applied in the 
form of a hat function over time. In particular, the traction vectors are given by 


256/9 t for Os <t< 2.5s, 
ta =—ty = f(t) | 512/9 | Pa, /f(t-—45-—t for 25s<t<5s, (493) 
768/9 0 for t> 5s. 


Table 4.3 provides a summary of the data used in the simulations. No Dirichlet boundary 
conditions are applied. 
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Material parameters H 997.5 Pa Geometry 
À 5209 Pa 
Specific heat capacity c 100 cha-* 
Expansion coefficient B 2283.10 K^ 
Thermal conductivity k 10 WK im^! 
Ref. temperature bo 293.15 K 
Mass density p 100 kgm~? 
Initial temperature da 290 K h 
0, 350 K 
Geometry h 10 m 
b 3 m 
Newton tolerance € 1075 - 
Simulation duration T 250 s 
Time steps At 0.08 s 
0.2 s 
0.4 s 


Table 4.3.: L-shaped block: Data used in the simulations 


Since in the initial loading phase the external forces are equilibrated, the total linear 
momentum of the block is a conserved quantity. In addition to that, after the loading 
phase (t > 5s) no external torque is acting on the block. 


071 104 


— L; (EM),-At = 0.4 — J, (EM),-At = 0.4 
— L, (EM),-At = 0.4 — J, (EM),-At = 0.4 
—— L; (EM),,-At = 0.4 — J, (EM),,-At = 0.4 


Figure 4.6.: L-shaped block: Total discrete linear momentum (EM), scheme (left), Total discrete angular 
momentum (EM). scheme (right) 
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Consequently, the total angular momentum is conserved as well. All of the integrators 
under consideration are capable to exactly conserve both momentum maps (up to numeri- 
cal round-off), independent of the chosen time step. This can be observed from Fig. 4.6, 
where representative numerical results of the energy-momentum integrator (EM), are 
shown. In the simulations three alternative constant time steps At € {0.08, 0.2, 0.4}s 
have been used. After the loading phase the total energy must be a conserved quantity. As 
expected, the (EM). scheme does exactly reproduce this conservation law (up to numerical 
round-off), see Fig. 4.7. In contrast to that, the schemes (M)g and (ME), are not capable of 
conserving the total energy. 


-105 
o [c 
3.35 
= 
v 3.3 
3.25 
0 50 100 150 200 250 


t [s] 
—— (EM),,-At = 0.08 
—— (EM),,-At = 0.2 


—— (EM),,-At = 0.08 
— (EM),,-At = 0.2 


— (EM),-At = 0.4 


—— (EM),-At = 0.4 


Figure 4.7.: L-shaped block: Total energy (EM). scheme (left), Incremental change of total energy (EM), scheme 
(right) 


It can be observed from Fig. 4.8 that both schemes have a tendency to increase the energy, 
depending on the time step. Typically, such energy blow-ups eventually lead to a failure 
of the iterative (Newton-Raphson) solution procedure. In the diagrams the break down of 
the simulation is indicated by vertical lines. Despite the capability of the (EM), scheme to 
conserve the total energy, the simulation typically breaks down at about the same point in 
time as the break down of the other two schemes occurs. The numerical instability of the 
(EM). is accompanied by a nonphysical decrease of the total entropy as can be observed 
from Fig. 4.9. In fact, the total entropy ought to be a non-decreasing function of time. Only 
the (ME), scheme is capable to correctly adhere to the second law of thermodynamics, 
independent of the time step (Fig. 4.10). Indeed, in each time step the total entropy does 
increase, as can be observed from Fig. 4.10 (right). The above investigations indicate that 
all of the three schemes under consideration are prone to numerical instability. 
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-10° 
3.5 2 
= 3.4 il = 
w | w 
3.3 | 
0 100 200 0 50 100 150 200 250 
t [s] 
—— (M),-At = 0.08 — (ME), -At = 0.08 
— (M)g-At = 0.2 — (ME),-At = 0.2 
— (M)g-At = 0.4 — (ME),-At = 0.4 
Figure 4.8.: L-shaped block: Total energy (M)g scheme (left), Total energy (ME); scheme (right) 
g g 
c c 
v v 


— (M),-At = 0.08 


— (EM), -At = 0.08 
— (M),-At = 0.2 


— (EM)„-At = 0.2 


— (EM)„-At = 0.4 — (M),-At = 0.4 


Figure 4.9.: L-shaped block: Total entropy (EM). scheme (left), Total entropy (M)g scheme (right) 
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— (ME), -At = 0.08 


— (ME), -At = 0.08 
— (ME), -At = 0.2 
— (ME),-At = 0.4 


— (ME),-At = 0.2 
— (ME),-At = 0.4 


Figure 4.10.: L-shaped block: Total entropy (ME),, scheme (left), Incremental change of total entropy (ME); 
scheme (right) 


After the loading phase, the environment of the present example can be characterized as 
thermally perfect in the sense of [146] (cf. Appendix A.2). Thus £ defined in (4.92) plays 
the role of a Lyapunov function that has to decrease with time. However, none of the 
schemes under investigation does correctly reproduce this behavior, as can be seen from 
Figs. 4.11 and 4.12. That is, depending on the time step and the duration of the simulation, 
all of the schemes inevitably exhibit numerical instabilities characterized by increasing 
values of £. 


—— (EM),,-At = 0.08 — (M)g-At = 0.08 
— (EM), -At = 0.2 — (M)g-At = 0.2 
— (EM),-At = 0.4 — (M)g-At = 0.4 


Figure 4.11.: L-shaped block: Lyapunov function (EM). scheme (left), Lyapunov function (M)g scheme (right) 
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-104 


100 120 140 160 180 200 220 240 
t [s] 
—— (MF)„-At = 0.08 
—— (ME),-At = 0.2 
—— (ME),-At = 0.4 


Figure 4.12.: L-shaped block: Lyapunov function (ME); scheme 


Eventually, the motion of the L-shaped block is illustrated in Figs. 4.13 and 4.14 with snap- 
shots at successive points in time. In addition to that, the distribution of the temperature 


and det (F) over the block is shown. 


Figure 4.13.: L-shaped block: Snapshots of the motion along with the temperature distribution over the block at 
t € (0, 32, 64, 96, 128, 160, 192, 224}s, obtained with the (M)g scheme and time step At = 0.08s 
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1.37 


0.95 


Figure 4.14.: L-shaped block: Snapshots of the motion along with the distribution of the determinant of the 


deformation gradient over the block at t € (0, 32, 64, 96, 128, 160, 192, 224] s, obtained with the (M)g scheme 
and time step At — 0.08s 


4.5.2. Rotating disc 


The second example deals with a rotating disc subjected to prescribed heat flow over part 
of the boundary surface (Fig. 4.15). The spatial discretization of the disc is based on 200 
tri-linear finite elements leading to a total of 360 nodes. 


Wo 1 
9 0.8 


0246 810 
t [s] 


Figure 4.15.: Rotating disc: Initial configuration and thermal boundary conditions (left), function f(t) for the 
prescribed heat flow over part of the boundary surface (right) 


The initial velocity distribution over the disc results from a prescribed angular velocity 
wo € R? and is given by 
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The initial temperature of the disc is homogeneously distributed and equal to the reference 
temperature 0p. In an initial period of time, t € [0, 4]s, heat flow is prescribed over one 
quarter of the lateral boundary surface (Fig. 4.15). In particular, the heat flow into the disc 
is described by 


sinit) for 0<t<4s, 


0 for t>4s. 


A plot of function f(t) can be found in Fig. 4.15. The rest of the boundary surface of the 
disc is assumed to be thermally insulated (¢ = 0). Note that the prescribed heat flow 
vanishes after t = 4s. For t > 4s, the complete boundary surface of the disc is assumed 
to be thermally insulated (q = 0 on 0,B = OB). Then, the environment of the disc is 
thermally perfect in the sense of [146] (cf. Appendix A.2). A summary of the data used in 
the simulation of the rotating disc can be found in Table 4.4. 


Material parameters À 3000 Pa Geometry 
u 750 Pa 

Specific heat capacity c 150 JK !m-^? 

Expansion coefficient B 1-1074 Kcti 

Thermal conductivity k 


20 WK^m^! 
Ref. temperature do 300 K <> "m 


Mass density p 893 kgm ? 
Radius r 08 m 

R 2 m 

t m 

€ 

T 


Thickness 0.4 
Newton tolerance 1078 - 
Simulation time 30 s 
Time steps At 0.04 s 
0.08 s 
S 


0.1 


Table 4.4.: Rotating disc: Data used in the simulations 
In the simulations constant time steps of size At € (0.04, 0.08, 0.1}s are applied. Since 
neither external loads act on the disc, nor displacement boundary conditions are imposed, 


the mechanical system at hand has translational and rotational symmetry. Consequently, 
the corresponding momentum maps are first integrals of the motion. 
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— L (ME), -At = 0.1 
— L, (ME), -At = 0.1 
— L (ME), -At = 0.1 


—J, (ME),-At = 0.1 


— J, (ME), -At = 0.1 
— J; (ME), -At = 0.1 


Figure 4.16.: Rotating disc: Total linear momentum (ME); scheme (left), Total angular momentum (ME); 
scheme (right) 


All three integrators at hand are capable to conserve the respective momentum map. 
Representative numerical results are shown in Fig. 4.16. 


-103 -1078 


E Wm n 1 


—— (EM),,-At = 0.04 — (EM),-At = 0.04 
— (EM), -At = 0.08 — (EM),,-At = 0.08 
—— (EM),-At = 0.1 — (EM),-At = 0.1 


Figure 4.17.: Rotating disc: Total energy (EM), scheme (left), Discrete change of total energy (EM), scheme 
(right) 
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-103 


am) - 
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— (M)g-At = 0.04 
— (M)g-At = 0.08 


— (ME),- At = 0.04 
— (ME), -At = 0.08 


— (ME),-At = 0.1 


Figure 4.18.: Rotating disc: Total energy (M)g scheme (left), Total energy (ME); scheme (right) 


Since heat flow into the system is prescribed in the initial time period [0, 4]s, the total 
energy is expected to increase. For t > 4s, the system is closed and the total energy should 
stay constant. Again the (EM), scheme is capable to correctly reproduce the first law 
(Fig. 4.17). However, despite the algorithmic energy conservation (for t > 4s), the (EM), 
scheme gets unstable, depending on the time step. The corresponding point in time of the 
break down of the simulation is indicated with a vertical line in the diagrams. At about 
the same points in time the other two schemes break down as well (Fig. 4.18). For these 
schemes the break down is accompanied by a sudden increase of the total energy leading 
to the divergence of the Newton-Raphson iterations. For the considered duration of the 
simulation (t € [0, 30]s), a time step of At = 0.04s is small enough to retain numerical 
stability of the three schemes at hand. 
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Figure 4.19.: Rotating disc: Total entropy (EM), scheme (left), Total entropy (M)g scheme (right) 
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— (ME), - At = 0.04 
— (ME),,-At = 0.08 


— (ME),-At = 0.04 
— (ME),- At = 0.08 
— (ME),-At = 0.1 


— (ME),-At = 0.1 


Figure 4.20.: Rotating disc: Total entropy (ME); scheme (left), Discrete change of total entropy (ME); scheme 
(right) 


Due to the prescribed heat flow into the disc, the total entropy of the disc is expected to 
increase in the initial time period [0, 4]s. For t > 4s, the system is closed and the total 
entropy should be a non-decreasing function of time. As expected, the (ME),, scheme 
is capable to exactly satisfy the second law, independent of the time step. This can be 
observed from Fig. 4.20. In particular, Fig. 4.20 (right) confirms that the change per time 
step of the total entropy is always positive. The (EM),, scheme does not correctly reproduce 
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the second law as can be seen from Fig. 4.19. Accordingly, the divergence of the iterative 
solution procedure is accompanied by an nonphysical decrease of the total entropy. The 
(M)e closely adheres to the second law as can be observed from Fig. 4.19. 


To shed further light on the numerical stability of the present schemes, we consider the 
Lyapunov function £ defined in (4.92). For t > 4s the system is closed and the function 
L£ should decrease with time (see also Appendix A.2). All three schemes fail to correctly 
reproduce this behavior, depending on the size of the time step and the duration of 
the simulation (Figs. 4.21 and 4.22). In particular, it can be seen that the smallest time 
step At — 0.04s yields a stable numerical simulation, at least in the considered time 
interval [0, 30]s. However, for larger time steps At = 0.08s and At = 0.1s, the numerical 
instability of each scheme becomes visible through the increase of £. 
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— (M)g-At = 0.04 
— (M)g-At = 0.08 


—— (EM),-At = 0.04 
—— (EM),,-At = 0.08 


— (EM),-At = 0.1 — (M),-At = 0.1 


Figure 4.21.: Rotating disc: Lyapunov function (EM), scheme (left), Lyapunov function (M)g scheme (right) 
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— (ME), - At = 0.04 
— (ME),-At = 0.08 
— (ME),-At — 0.1 


Figure 4.22.: Rotating disc: Lyapunov function (ME), scheme 


In addition we investigate the order of accuracy of the present schemes. We therefore intro- 
duce the relative error in the position, linear momentum density and absolute temperature 


defined by 


tol 


.. I-P. _ IIP- P. Ilz: — 10—9.llz = 
eo pl, ^ ep Plz, > eo Ma lle |l. = [fs (e. 0) aV] ; 


where index r refers to the reference solution obtained with a very small time step. 
Moreover, (e, e) denotes the inner product. As can be observed from Figs. 4.23 and 4.24, all 
schemes are second order accurate in the position, linear momentum density and absolute 


temperature. 
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Figure 4.23.: Rotating disc: Error in the position (left), Error in the linear momentum density (right) 


10-9 


Figure 4.24.: Rotating disc: Error in the absolute temperature 


Eventually, the motion of the disc is illustrated in Fig. 4.25 with snapshots at successive 
points in time. In addition to that, the distribution of the temperature over the disc is 
shown. 
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Figure 4.25.: Rotating disc: Snapshots of the motion at successive points in time t € (0, 4, 8, 12, 16 ,18 ‚24, 
28}s, and corresponding temperature distribution, calculated with the (M)g scheme and At = 0.04s 


4.5.3. Rotating disc in a thermally perfect environment 


The last example of this Chapter deals with a slight modification of the previous example 
which renders a thermally perfect environment for the disc throughout the whole duration 
of the simulation (Fig. 4.26). 
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Figure 4.26.: Rotating disc in a thermally perfect environment: Initial configuration, temperature boundary 
conditions imposed on a quarter of the lateral surface of the disc and initial distribution of the temperature 


Instead of prescribing the heat flow over a quarter of the lateral surface of the disc, we 
now impose temperature boundary conditions on the same portion of the lateral surface 
of the disc. In particular, the temperature at the finite element nodes located on the 
boundary in question is constrained to be equal to the reference temperature 99. The 
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initial temperature at the remaining nodes of the finite element mesh is chosen to assume 
the value du + 80K. The remaining simulation data are identical to those used in the 
previous example (see also Table 4.4). A summary of the data used in the simulation of the 
rotating disc in a thermally perfect environment can be found in Table 4.5. In this example 
we focus on the temperature-based (M)g scheme, which makes possible to impose the 
temperature boundary conditions in a straightforward way through standard Dirichlet 
boundary conditions. Constant time steps of size At € (0.04, 0.075, 0.0875 }s are applied. 
Again the mechanical system at hand has both translational and rotational symmetry. 


Material parameters 


Specific heat capacity 
Expansion coefficient 
Thermal conductivity 
Ref. temperature 
Mass density 

Radius 


Thickness 
Newton tolerance 
Simulation time 
Time steps 


3000 
750 

150 
1.1074 
20 

300 
8.93 

0.8 

2 

0.4 
10-8 
70 

0.04 
0.075 
0.0875 


Pa Geometry 


5 “> T 


nnnn 


Table 4.5.: Rotating disc in a thermally perfect environment: Data used in the simulations 


The corresponding momentum maps are exactly conserved by the (M)g scheme as can be 


observed from Fig. 4.27. 
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Figure 4.27.: Rotating disc in a thermally perfect environment: Total linear momentum (M)g scheme (left), Total 


angular momentum (M)g scheme (right) 


Due to the temperature boundary condition heat is flowing out of the disc. The corre- 
sponding evolution of the total energy is depicted in Fig. 4.28 (left). In addition to that, 
the evolution of the total entropy is shown in Fig. 4.28 (right). 
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Figure 4.28.: Rotating disc in a thermally perfect environment: Total energy (left), Total entropy (right) 
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Figure 4.29.: Rotating disc in a thermally perfect environment: Lyapunov function 


As before, the break down of a simulation due to numerical instability is indicated by 
vertical lines. Again the onset of numerical instability can be detected from the plot 
of the Lyapunov function in Fig. 4.29, although this time the nonphysical increase is 
not as pronounced as in the two previous examples. Eventually, both the motion and 
the evolution of the temperature distribution are shown in Fig. 4.30 with snapshots at 
successive points in time. 


380.00 


353.33 


326.67 


300.00 


Figure 4.30.: Rotating disc in a thermally perfect environment: Snapshots of the motion at t € (0, 10, 20, 30, 
40, 50, 60, 70}s and corresponding temperature distribution, calculated with the (M)g scheme and At = 0.04s 
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5. Energy-momentum-entropy 
consistent numerical methods for 
large-strain thermo-elasticity 
based on the GENERIC formalism! 


The results of Chapter 4 lead to the conclusion that all major balance laws should be 
maintained in the discrete setting in order to enhance the numerical stability of numerical 
methods for finite-strain thermo-elastodynamics. Therefore fully structure-preserving 
numerical methods are developed in the present chapter. We will resort to the GENERIC- 
based description of large-strain thermo-elasticity, which makes the free choice of the 
thermodynamic state variable possible. In particular, one may choose (i) the internal energy 
density, (ii) the entropy density, or (iii) the absolute temperature as the thermodynamic 
state variable. Three alternative energy-momentum-entropy (EME) schemes result from 
the present approach that make use of a GENERIC-consistent space discretization. These 
schemes are directly linked to the respective choice of the thermodynamic state variable. 
Numerical examples confirm the structure-preserving properties of the newly developed 
(EME) schemes, which exhibit superior numerical stability. 


5.1. Large-strain thermoelasticty 


In this section we briefly recapitulate the variational formulation of large-strain thermo- 
elasticity with heat conduction of Chapter 4 which lies at the heart of the proposed 
discretization in space and time. In contrast to Chapter 4 we introduce a frame-invariant 
formulation from the outset. 


5.1.1. Underlying variational formulation 


We consider a continuum body with material points X = X ^e, in the reference configu- 
ration B C IR?, see Fig. 5.1. Here and in the sequel the summation convention applies to 
repeated indices. Moreover, e4 denote the canonical base vectors in R3. 


1 This Chapter is based on [2] 
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€i 


€3 


Figure 5.1.: Reference configuration B with boundary OB and current configuration (B, t) at time t. External 
tractions t = FSN act on the boundary of the current configuration. In addition to that, the heat flux across 
the current boundary is denoted by q = Q - N. 


Within the Lagrangian description of continuum mechanics the deformed configuration 
of the body at time t is characterised by the deformation map : B x T +> R°, where 
T = [0, T] is the time interal of interest. The velocity of the material point X € B located 
at x = (X, t) is given by v = ¢, where a superposed dotes denotes the material time 
derivate. The deformation gradient is given by 


F = ôx. (5.1) 


Main ingredients of the GENERIC framework are the internal energy and the entropy. In 
addition to that, the choice of the thermodynamic state variable plays an important role. 
Following Chapter 4 we allow for the free choice of the thermodynamic state variable 
T : B x I + R from among three options 7 € (0,7, u). These options are: (i) the 
absolute temperature 0, (ii) the entropy density 7, and (iii) the internal energy density u. 
Depending on the choice of the thermodynamic state variable, the absolute temperature 
can be written in the form (see also [94] and [100]) 
|. Ou (C, 7) 


e'(C,7) = 37 (6,7) . (5.2) 


Now the internal energy density and the entropy density, respectively, are given by 


u =u (C,T) and qnm. (5.3) 
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In this connection, a frame-indifferent constitutive formulation for thermoelastic materials 
is based on the right Cauchy-Green tensor C — FTF. The GENERIC-based weak form 
pertaining to large-strain thermo-elasticity can be written in the form (see Chapter 4) 


o-f we: (p-p) dV, 
o=) (orp (BED Bien | wasted, 
B 0,8 
R =]; Wr / Wr Wr _ 
= 2 : — —— |: d dA 
0 f (sevo p) (a5 ) v(zz) a) ed 
(5.4) 
where p : B +> R+ is the mass density in the reference configuration. Moreover, p : 
B x I++ R? is the linear momentum density and b : B +> R? represent prescribed body 


forces which, for simplicity, are assumed to be dead loads. The second Piola-Kirchhoff 
stress tensor is given by 


S= S'(C, T) L2 (cu = 0’0cn’) : (5.5) 


Furthermore, the material heat flux vector assumes the form 


Q = Q(C,7) = (6?K'v (&) (5.6) 
where K = K'(C, 7) is the positive semi-definite material conductivity tensor. The weak 
form needs be supplemented with initial and boundary conditions, respectively. For that 
purpose, the boundary OB of the continuum body is decomposed into a displacement 
boundary 0,8, on which » = Ẹ is prescribed, and a traction boundary 0,B, on which 
the external traction t is prescribed such that FSN = t (Fig. 5.2). In this connection, the 
standard relations 0,6 U 0,B = OB and 0,B N 0, B = () hold. Similarly, for the thermal 
part we consider the subsets 0,6 and 0,B, with the properties d,B U ôB = OB and 
-BN 0,B = 0 (Fig. 5.3). The thermodynamic state variable is prescribed on 0,8, i.e. 
T = T, whereas the heat flux is prescribed on 0,B, i.e. Q- N =q. 


The unknown fields are subject to initial conditions of the form y(-,0) = X, p(-,0) = pvo, 
and 7(-,0) = m in B. Here, vo is a prescribed velocity field and ro is a prescribed field of 
the thermodynamic state variable r € (0,5, u}. 
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€3 


Figure 5.2.: Mechanical part of the IBVP. Note that t = FSN denotes prescribed external Piola tractions acting 
on the current boundary expressed per unit area of the reference boundary 05 B. 


ea 


Figure 5.3.: Thermal part of the IBVP. Note that = Q - N is the prescribed rate of heat transfer accross the 
current boundary expressed per unit area of the reference boundary 948. 


5.1.2. Frame-invariant local form of the field equations 
We next deal with the frame-invariant local form of the field equations corresponding 


to the GENERIC-based weak form (5.4). In this context, we consider alternative frame- 
invariant formulations related to the specific choice of the thermodynamic state variable 
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T € (0,n, u}. Using standard arguments from the calculus of variations, weak form (5.4) 
gives rise to the corresponding local form of the field equations 


REESE 
gu (5.7) 
p = b + Div (FS) , 
together with 
EE di nt Di 
T= B cu :V(p pP) - Bul ivQ. (5.8) 


These equations have to hold for all X € B and t > 0. The last equation governs the 
time-evolution of the thermodynamic state variable r € (0,1, u}. Note that the second 
Piola-Kirchhoff stress tensor S in (5.7)? assumes the form (5.5), and the material heat flux 
vector Q) in (5.8) is given by (5.6). 


Choosing T = 0 leads to the temperature-based formulation, where the second Piola- 
Kirchhoff stress tensor follows from (5.5) and assumes the form 


S = 2[0cu(C,0) — ðc (C, 0)] . (5.9) 
The temperature typically enters the definition of the Helmholtz free energy density 
V — v(C,0). (5.10) 
Then the entropy density is defined by 
n =7(C, 0) = -0jJ (C, 0) . (5.11) 


The internal energy density function then follows from a Legende transformation of the 
free energy density. Accordingly, the internal energy density is given by 


u — u(C,0) = v(C,0) + 05(C,0). (5.12) 


Taking into account the temperature formula (5.2) along with Y% = u — 67), it can be seen 
that expression (5.9) for the second Piola-Kirchhoff stress tensor is equivalent to 


S = S(C,0) = 20cU(C,0) . (5.13) 


Equations (5.11) and (5.13) are sometimes called thermal equations of state (see, for exam- 
ple, [147]). 


Alternatively, we may use the entropy density 7 : B x Z ++ Ras thermodynamic state 
variable leading to the entropy-based formulation. Accordingly, choosing 7 = n, the second 
Piola-Kirchhoff stress tensor follows from (5.5) and assumes the form 


S = S(C, n) = 28cu(C, n) . (5.14) 
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free energy internal energy entropy temperature 
T v'(C. 7) w (C. T) i! (C, 7) 9'(C.7) 
0 ¥(C, 6) u(C, 0) n(C, 8) 
n v(C, n) u(C, n) n e(C. n) 
u Y(C, u) u iC, u) O(C, u) 


Table 5.1.: Notation used in the present Chapter, depending on the choice of the thermodynamic state variable 
T € {6,n, u}. 


Temperature formula (5.2) leads to 
O(C, 7) = (C, n) . (5.15) 


The last two equations belong to the caloric equations of state (see, for example, [136]). 


Yet another option is to use the internal energy density u : B x Z ++ Ras thermodynamic 
state variable leading to the internal-energy-based formulation. Thus, choosing 7 — u, the 
second Piola-Kirchhoff stress tensor follows from (5.5) and assumes the form 


S = $(C,u) = —20067(C, u) . (5.16) 
Furthermore, temperature formula (5.2) yields 


^ 1 

O(C,u) = — =. (5.17) 
9f (C, u) 

The notation used in the present work is summarized in Table 5.1. A more detailed 

treatment of the alternative choice for the thermodynamic state variables can be found in 

[100, Sec. 2.3] (see also Section 4.2). 


5.1.3. Balance laws 


We next consider important balance laws of the underlying continuous formulation of 
finite-strain thermo-elasticity. The main goal of the present work is to preserve these 
balance laws under discretization. In this section we focus on the pure Neumann problem 


(Le. ðB = 0,B = OB). 


The linear momentum of the continuum body is defined by L = Js p dV. Inserting 
Wp = € into (5.4)5, where & € IR? is arbitrary but constant, leads to 


dL 
£y —— ( [pav - [ rswaa) (5.18) 
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Due to the arbitrariness of £ € R?, the last equation recovers the balance law for linear 
momentum. The parentheses on the right-hand side of (5.18) contain the resultant external 
load applied on the continuum body (see also Fig. 5.1). 


The total angular momentum relative to the origin of the inertial frame is defined by 
J= jJ q x pdV. Substituting w, = p x 6, wp = E x v, into (5.4) and subsequently 
adding the equations yields 


dJ 
em e (fosoave [ oxrsnad) . (5.19) 


Due to the arbitrariness of € € R®, the last equation recovers the balance of angular 
momentum. Note that the parentheses on the right-hand side of (5.19) contain the resultant 
external torque about the origin (see also Fig. 5.1). 


For verification of the balance of energy we choose the following test functions in (5.4) 


= ! E 
Wy = p€ , wo=-b, 
ge — „1 
Wp = p€ , or Wp-—p P, 
1 
Wr = ore , Wr — Oru’ , 


where the total energy density e’ is given by 


1 
(p, F, pT) = 5p 'p:ptu(C,r) b.e. (5.20) 
Substituting these quantities into (5.4) and subsequently adding the three resulting equa- 


tions, a straightforward calculation recovers the balance law for energy 


dé’ 


= a. um 5 
anh. (po 'p-FSN-—N-Q) dA. (5.21) 


Concerning the balance of entropy, we insert w, = 0,7 into (5.4)3, to obtain 


o-f 0,n't + V(p !p): 2F8cy) - V ar -Q av+ f Ont Q- NdA 
B d P : enu TA ap Oru’ 


ffl a LP 
- [ (35 -e (5) oie (2) av f banaa. 


Here, use has been made of formula (5.2) for the temperature along with expression (5.6) 
for the material heat flux vector. Moreover, the identity ¢ = p^ ! p has again been taken 
into account. The above equation can be rewritten as 


ds! 1 > 1 1 

—- — |. K d -NdA. .22 

di Ir (s) i v(s) T a P 
Sa a_r 


>0 


The last equation corresponds to the Clausius-Duhem form of the second law of thermo- 
dynamics (see, for example, [136, Sec. 5]). 
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5.2. Discretization in space 


Next, we perform the discretization in space of the GENERIC-based weak form (5.4). To 
this end we apply the isoparametric finite element approach (see, for example, [145]), 
based on finite-dimensional approximations of the following quantities 


p"(X,t) = N"(X)q,(t,  v°(X, t) = N*(X) va(t), (5.23) 

and 
Th(X,1) = N*(X) ta(t). (5.24) 
Here, the summation convention applies, where a = 1,..., N, and N denotes the total 


number of nodes in the finite element mesh. Moreover, N^ : B — Rare the nodal shape 
functions and q, (t), v,(t) € R°, Ta(t) € R are the respective nodal values at time t. 
Analogous approximations are used for the test functions w,, Wp and w,, denoted by 
wh wp and wP. Then weak form (5.4) leads to the following semi-discrete equations: 


p? 
o= fwg (^ - v") dV, 
= f (ws: (on - b) evene) av- | wb-taa, 
9 7 (5.25) 
h -h h h h 
= | wh [55 veh: | — Fha dV 
i: ( (Mo en) 


- fv acis ) grav + f UNE. od 
"RINCNPY 9,8 Hn(O ub) 4 
where 


sh =2 (Ocu® eo den‘) ; 


: 1 
Q" = (Kv (ax) i (5.26) 
h_ II; (u^) 
I‘ 


In this connection, ut = w(C^, 7h), jh = w(C^, 7^), and K" = K'(C^, 7^). The 
interpolation formulas in (5.23) give rise to 


F'—q &VN* and C"—-q,.q, VN" & VN*. (5.27) 


Moreover, II, (0,u") denotes the Lə projection of function O, u^ into the finite-dimen- 
sional space spanned by the shape functions N^, a — 1,..., N. That is, 


II; (0,u") = N*(8-u)a , (5.28) 
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where the nodal values (0,u)q are determined by 
1 N“ (8.u" — II, (3-u*)) av =0, (5.29) 
B 


for a = 1,..., N. In particular, (5.29) together with (5.28) constitute a linear system of 
algebraic equations given by 


H^ (0.u), = ] vo dV , (5.30) 
B 


where H®? denote the components of the positive definite Gram matrix [H^] defined 


by 
H® = if N°N? av. (5.31) 
B 
We note for later use that (5.28) in conjunction with (5.30) lead to the relationship 
IIj(O.uP) = N° Ha; / N*8.u^ dV , (5.32) 
B 


where Ha, denote the components of the inverse of the Gram matrix, ie. [Hay] = 
[H= 


Analogous relationships hold for the projection of function 0,7>, IT, (8-n™). It is worth 
mentioning that the projections IT),(0,u") and Hp (0, r^) are only required if the func- 
tions 0,u> and 0,7" do not belong to the finite element space spanned by the shape 
functions N“. For example, in the temperature-based formulation (ie. for 7 = 6), 
dou® = dou(C", 0^) corresponds to the specific heat at constant deformation. Thus, if 
this quantity is prescribed to be constant, the projection Il, (09%) need not be performed. 
However, in general the temperature-based formulation requires both projections (i.e. 


II; (Og) and II; (8o7])). 


In contrast to that, the two alternative formulations based on the choice 7 € {n,u} 
in general require only one projection. In particular, the formulation in terms of the 
internal energy density relies on II, (0,17), whereas the entropy-based formulation relies 
on II; (ù). Originally, the projection has been introduced in the framework of the 
entropy-based formulation in [102] (see also [106]). 


In the next Section we will show that the proposed discretization in space of the GENERIC- 
based weak form (5.4) yields semi-discrete equations which can be brought into GENERIC 
form for finite dimensional systems (2.1). Thus, we speak of a GENERIC-consistent space 
discretization, see Section 2.1.2. In this connection the introduction of the above projections 
is essential for retaining consistency of the semi-discrete formulation. 
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5.2.1. GENERIC-consistent space discretization 


In the following we will show that the proposed discretization in space of the GENERIC- 
based weak form (5.4) yields semi-discrete equations which can be brought into GENERIC 
form (2.1), thus we show that the above proposed space discretization is GENERIC- 
consistent in the sense of Section 2.1.2. 


The present interpolation of the velocity field gives rise to nodal velocity vectors vq, 
cf. (5.23). We first introduce the conjugate nodal momentum vectors via a standard 
Legendre transformation. Accordingly, we express the kinetic energy Ekin = &,(p) = 
Js ip |p : pdV in terms of the velocity field by using the identity p = pv, leading to 
Exin = T(v) = fniov -v dV. Using the discrete velocity field v^ yields the kinetic 
energy of the discrete system given by 


1 
PETS gM va Vp, (5.33) 


where the components of the standard consistent mass matrix [M""] are given by 
Me? = / pN“N?aV. (5.34) 
B 
Note that the mass matrix [M""] is symmetric and positive definite. Now the conjugate 
nodal momentum vectors are defined by 
p° = ôv T = M**y, . (5.35) 


The last equation implies that the nodal velocity vectors can be written in terms of the 
nodal momentum vectors via 


Va = Map’, (5.36) 


where Ma» are the components of the inverse mass matrix, i.e. [May] = [M?*]-!. Now 
we introduce the state vector of the discrete system at hand which is comprised of the 
nodal quantities q,, p^, Ta, a = 1,..., N. We arrange these quantities in the nodal state 
vector of the system given by 


z = (qı, 605 qy; pt, e Poo T +. TN) * (5.37) 


Hereby the state vector has to be viewed as a column vector. Next we aim at the total 
energy and the total entropy of the discrete system which play the role of generators in 
the GENERIC formulation (2.1). The two required functions assume the form 


1 
€ (z) = [wav+ ; Map" p'-q,: | N°bav , 
B B 


sta) = [ av. 


(5.38) 
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Concerning the derivatives 0,€ and 045, we obtain the individual contributions 


ða E = [vn av - f Nebav, 
B B 
pE = Mapp’ , (5.39) 
0 = nor dv , 
B 
and 
ðq S = / Ogg" VN" dV , 
B 
Op. S — O0, (5.40) 


0, S = ] von dV . 
B 


Next, we take a closer look at the semi-discrete formulation emanating from the space- 
discrete weak form (5.25). To this end we focus on the pure Neumann problem and neglect 
the boundary integrals in (5.25). This is in accordance with the fact that the basic form 
(2.1) of the GENERIC formulation is valid for closed systems. The kinematic equation 
(5.25) directly leads to H?^(q,, — Va) = 0, where the components of the Gram matrix 
[H^] have been introduced in (5.31). Due to the positive definiteness of the Gram matrix 
one obtains 


Gq = Mapp? = Opa. (5.41) 
The space-discrete weak form of the balance law for linear momentum, (5.25)5, yields 
Mey, = / N“bdV - / F's'vwN?qy . (5.42) 
B B 
Taking into account (5.35), (5.26)1,; and 


ru? = 2F cul, 


den = OF "den . (5.43) 
into account eq. (5.42) can be recast in the form 
p° = i N*bdV — / Ogu^ V N? dV + T Ta (ru). n Na dv. (5.44) 
B B gll (Or m") 
With regard to (5.32) and (5.39)3, we get 
TI,(0,u>) = NH opôr E . (5.45) 
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Accordingly, taking into account the last equation and (5.39),, (5.44) can be rewritten as 


N4 
p^ — IE + / Tor v^ dV HajO,, E . (5.46) 
B T: 


The space-discrete weak form of the energy equation, (5.25)5, yields 


Hg = — Vd y ay h Ogg? VN7 dV 
g Un(0;7") (5.47) 
LV N°? (9")?K'V II, (0. ) dV 
II, (ð uh) II, (ð uh) 


In analogy to (5.45), we have 
IIj(8.) = N*Ha0,6. 
Using the last equation together with (5.36) and (5.39)5, (5.47) can be recast in the form 


N® 
ta = — HavOpaE - | ————~ Ir V NI 
e moa” dd 


+H fy PE M (G^? K^v M L| dVHaO,,8 
ab 5 II; (0, u^) II, (8, uh) cdOr49 - 


To summarize, the semi-discrete evolution equations resulting from the discretization in 
space of the underlying GENERIC-based weak form can be written as 


(5.48) 


a = pa E ? 
p^ = —Og,€ + 1°, 0r,€ , (5.49) 
PAT 
ee (v) 8p£ + MarInS , 
where 
N’ 
1%; = hi UN V Ha, 
B in AOT) (5.50) 


— H, fy MU a (Gy? K^v TN dV H, 
Mab = Iac a II; (9, u^) II; (9, u^) db + 


The evolution equations (5.49) fit into the GENERIC framework for finite-dimensional 
systems provided by (2.1). In particular, (5.49) can be written in the form 


z=L0,E -MOS, 
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where the state vector is given by (5.37). Moreover, the Poisson matrix reads 


E 0 KA s (5.51) 


0 Eok 0 


where I is the identity matrix (with appropriate dimension corresponding to the partition- 
ing of the state vector (5.37)), and matrix [1^5] consists of vectors 1”, defined in (5.50)1. 
Specifically, we have 


la NA p 
LS e i [e (5.52) 
iNi k 
The friction matrix is given by 
0 0 0 
M= 0 0 0 ; (5.53) 
0 0 [mab] 


where the coefficients m.» have been introduced in (5.50)2. It can be easily verified by 
a straightforward calculation that the two degeneration conditions (2.2) and (2.3) are 
satisfied. Thus the validity of (2.4) and (2.5) ensures thermodynamic consistency of the 
semi-discrete formulation. 


Remark 6. Instead of using the interpolation of the velocity field (cf. (5.23)), one could as 
well interpolate the linear momentum density field. This would lead, for example, to minor 
changes in the partial derivative (5.39)5 and the Poisson matrix (5.51). However, it would 
not affect the thermodynamic consistency of the resulting semi-discrete formulation. 
These observations are valid for constant mass density p. We prefer the interpolation of 
the velocity field, because (i) it does not put any restrictions on the mass density and, (ii) it 
leads to nodal velocities as degrees of freedom on which initial and boundary conditions 
for the velocity field can be directly imposed. 
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5.2.2. Balance laws 


We next give a short summary of the balance laws pertaining to the space-discrete for- 
mulation described above. We restate the semi-discrete evolution equations derived in 
Section 5.2.1 in the form 


d = Va ’ 
Mew, = / N*bdV — / Ogu^ VN" dV + if larn VN* dV , 
B B B 
N^ N^ 
Hh = -v | Opn V N? dV - [v (mos) .Khvo^qv 
EA i cT g^ MI, (0r u*) 


N° 
gdA. 
Lam 


(5.54) 


For simplicity of exposition we have neglected in (5.54)2 the standard contribution of the 
external boundary tractions. The total linear momentum of the discrete system at hand is 


given by 
un = f ov*av 
B 
= f pN* dVv, 
B 
N anrb 
= jr N’ dVv, 
= JB 


5 N ab 
= M Vb. 
a=1 


Note that use has been made of the completeness condition RU N“ = 1 for the nodal 
shape functions. Now (5.54)? directly leads to the balance law for linear momentum 


(5.55) 


Tipa / bdV, (5.56) 


which can be regarded as the semi-discrete counterpart of (5.18). The total angular 
momentum relative to the origin of the inertial frame is given by 


a= f op x v” dV 
B 


(5.57) 
= M**q, X Vp. 
A straightforward calculation based on (5.54); 2 yields 
d n h 
&£.—J'—£. | e^ xbdv, (5.58) 
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for arbitrary but constant € € R°. In the last equation the following relationship has been 
taken into account: 


E- (da X (ðpu” — O^Ogy^) VN") ((Oru" — G^ügg")V N* Q qa) 
((aeu" = 6" pn )F*") 


é (Frst rt") 
0. 


| 


E: 
E: 


(5.59) 


Here, use has been made of (5.43) and (5.26),. Moreover, € is skew-symmetric such that 
£a = £x a for any a € R°. Relationship (5.58) can be viewed as semi-discrete counterpart 
of the balance law for angular momentum (5.19). 


For completeness of exposition we eventually verify the balance laws for energy and 
entropy, although the GENERIC-consistent space discretization at hand guarantees ther- 
modynamic consistency as has already been shown in Section 5.2.1. Multiplying (5.54)3 
by (0,u)q, and taking into account (5.28) together with (5.26)3 yields 


(0,u)gH 7 = -va J O^8gnp^v N? dV — jdA. 
B 04B 


With regard to (5.30), the last equation can be rewritten as 
f O,u ^ dV + f Oden: F” av = — I jdA. (5.60) 
B B 8,B 
Scalar multiplication of (5.54) by v, yields 
Va Mev, = v,- ( f N°bdV — / ru” VN" dV + / Oden VN" av) 
B B B 


or 


1d 
5g (ve Mw) = f 


v bdy = foe E^ av + [ eos Ev. 
B B B 


Substituting from the last equation into (5.60) leads to 


. -h ld 3 
[low topu: È ) dV + 2a uM vv) - f, 


The last equation can be rewritten as 


d 
—£gh- -f jdA, (5.61) 
dt 5,8 
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where the total energy of the discrete system assumes the form 
1 
E! = n dV 4-2M?hy,. vy — je -bdV. (5.62) 
B 2 B 
Note that (5.61) can be regarded as semi-discrete version of the balance law for energy 


(5.21). 


To recover the balance of entropy, multiply (5.54)3 by (0-7), and take into account (5.26)3, 
to obtain 


1 1 
(O-n) gH” iy = —va- ji Ogg" VN? dV — / vV(<):K?veorav - N — dA. 
B B oh 04B e 


The last equation can be rewritten as 


h-h h.p” 1 hy2qgh 1 l- 

- :F )dV = — |: K d dA, 
or + FN ) V lv x) (6?) v(a) V - ond 
SS ES 9€ — ( 

=45h >0 
(5.63) 
where the total entropy of the discrete system is defined by 
Sh- jr dV = faha. (5.64) 
B B 


Note that (5.63) can be viewed as semi-discrete version of the balance of entropy (5.22). 


5.3. Discretization in time 


We aim at a second-order accurate, implicit time-stepping scheme which is capable of cor- 
rectly reproducing the main balance laws outlined above. Due to its structure-preserving 
properties this type of integrator is called energy-momentum-entropy scheme. To devise 
such a scheme, we essentially apply the mid-point rule in which the derivatives of the 
internal energy density and the entropy density, respectively, are replaced by appropriate 
discrete derivatives. 


We focus on a representative time interval [t,, tn+1] with corresponding time-step size 
At = tn+1 — tn. The discrete approximations at times t» and tn+1 of a function f(t) will 
be denoted by fn and fn+1, respectively. Assume that the nodal state variables at time 
tns qa,» Vans and Ta, are given. The associated fields result from the nodal interpolation 
formulas (5.23) and are denoted by ql, v^ : Br R? and zb : BH R, 7 € (0,9, u}. 
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5.3. Discretization in time 


We aim at the determination of the corresponding quantities at time £541 through the 
mid-point type discretization of the semi-discrete formulation (5.25) given by 


0= wh. amp +Vwhs F? sh dV — w! -t,,1dA 
B P At p Emh algo 5,8 p n+5 , 
h h 
_ hf n+1 7 7n l^ x 2 h h 
o= [wt ( Apo Yen (morso )) 4 


fy 2 Qh av+ f Wr og dA 
" II; (Dur) algo ] 9,8 II; (Dub) +3 , 


(5.65) 
where 
Sigo —-2(Dcu^ — Oh. Den") ; 
1 
h | (gh 23h 
Qalgo = (8s) KV (a) 2 (5.66) 


oh = II; (Dut) 
algo II; (D;7") E 
h | «(e^ h 

and Kj,, = K (Coys Taga): 

The above discretization in time relies on the use of discrete derivatives in the sense of 
Gonzalez [46]. In particular, the discrete derivatives are applied to the internal energy 
density and the entropy density, respectively. For example, in the case of the internal 
energy density, the discrete derivatives are denoted by Dyu” and Dou, respectively. In 
particular, D,u" is defined by 


1 
Dru = 5 (dugs (à rta) dug. (Ar), (5.67) 
where 
wu (C, m41) — v (C, m) - O-u'(C, Tny) AT 


duc (Tn, Tn41) = ru (C, 7,41) + (Ar)? i 


and Ar = 7544 — Tn. Furthermore, Deu" is defined by 
1 
Dou" = 2 (du, (Cb, Cu) + du, n (Ch, Gea) , (5.68) 
n n+l 
where 
du+(Cn, Ca+1) = Ocul (Cangi T) 


u'(Cy41,7) — v (Cy, 7) — Icou'(C, 41,7) >: AC NG 
AC: AC i 
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and AC = C,,4; — Cy. It can be verified by a straightforward calculation that the discrete 
derivatives (5.67) and (5.68) satisfy the directionality condition 


n n? n 


Deut : (Chis — ci) + Dub Ua — TR) = WO qm —w(Ch, zh). (5.69) 


Analogous considerations apply to the discrete derivatives of the internal entropy density, 
D-n? and Dan”, respectively. Moreover, the time-average of any quantity (e) is given by 


z((*). 4 (nr): In particular, this implies 


(Ca Eg Cu) i (5.70) 


Note that in general C,, 1 # Foyt Far} In (5.65) and (5.66), II;,(D;u") denotes the L2 


projection of Du! into the finite dimensional space spanned by the shape functions N“, 
a = 1,..., N. In analogy to (5.30), 


H** (Du), = ] ro dV , (5.71) 
B 


constitutes a linear system of algebraic equations for the determination of the nodal values 


(Du), leading to 
TI,(D-u") = N*(Dzu), . (5.72) 
Analogous relations apply for II; (D... 7^). 


Remark 7. The standard mid-point rule can be recovered from (5.65) by simply replacing 
the discrete derivatives with the mid-point evaluation of the corresponding standard 
derivatives. That is, in the mid-point rule, instead of (5.67) and (5.68) one has to choose 


ho 9 (uhT ph h h_ IphT ph h 
D-u = O-w (Fr Far Tati) and Deu = deu (Fr Fati Taya) - 


5.3.1. Discrete balance laws 


We next show that the present discretization in space and time does indeed yield an EME 
scheme. In particular, the discrete fulfillment of the main balance laws can be shown along 
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5.3. Discretization in time 


the lines of the semi-discrete formulation dealt with in Section 5.2.2. To this end we recast 
the scheme (5.65) in the form 


da, m qda, _ 
At ath? 
ab Ybn+ı — Vbn a h a h h a 
M = | N*bav— | Deu" VN* dV + | 8, Dg" V.N* dV , 
At B B B 


ab Tbng1 — Ton _ A hy yd 


Ne 5 " Ne u 
-K d dA 
n (ow) algo V alee V T II; (Dun) +2 ’ 


(5.73) 


where Oligo has been introduced in (5.66)3. Moreover, similar to (5.43), 


Dru = oF , Deu? ; 
2 


(5.74) 
Den" = >) yt Der . 


Note that again the standard contribution of external tractions has been neglected in 
(5.73) for simplicity. Since the total linear momentum is given by L” = 51 My, it 
immediately follows from (5.73)5, that 


Lb -Lh-At / bdV. (5.75) 
B 

The last equation reflects the consistent approximation of the balance law for linear 

momentum. With regard to (5.57), the total angular momentum relative to the origin of 

the inertial frame assumes the form J” = M abq | X Vp. Then we have 


ees ee ad on X (Vs; — Von) + (Gang: — Aan) <a) 6.76) 


Substituting from (5.73), » into the last equation leads to 


E- (Ja Ja) - € (a... x M” (vor — vs,)) 
i (5.77) 


=E At | oh, xbav, 
B 


for arbitrary but constant € € R. In the last equation it has been taken into account 
that 


E. ce x (Deu^ — eh Der VN") 


| 


É: ((Dru" — Ot Des") FH,) 


P ES h ght (5.78) 
er (Pha Sie FT) 
0. 


| 
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This result follows from the symmetry of See along with the skew-symmetry of 3 


Note that (5.77) confirms the consistent approximation of the balance law for angular 
momentum. 


To recover the discrete version of the balance of energy, multiply (5.73)3 by (D-u). and 
take into account (5.72) together with (5.66)3 to obtain 


(Druja Alta — Ton) = —Atva Í Oh, Der" VN? dV — At I Tapi dA. 
nts B ə B 2 


q 


Employing (5.71) and (5.73)1, the last equation can be rewritten in the form 


J 5e — rb) qv + J 9. oe : (FR, — F?) dV = -At | Tny} dA. 
B B 0,8 


q 


(5.79) 


Scalar multiplication of (5.73)2 by Va M yields 
ntj 


ab Vbn+ı — Vb, — . an _ h a 
pr uM CUAL Na [xe Dgu" VN“ dV 


Va 


/ Oh Den VN* dV , 
2 B 


or 


Me 
or (Vangi © Vbagi — Van * Vb.) =At | vua -bdV — [pr : (FR, - Fav 


4 rer : (Fb, - Fav. 
B 
Substituting from the last equation into (5.79) and taking into account the relationship 
Dru" : (Far -F,) = 2F5, Deu : (Far =E) 
= Dou": (om = Ch) , 


yields 


: à Me 
I (D^ (ris A Th) + Deu" : (CER 73 ch) dV + 9. (Vasa “Vagi = Van € Von) 


= f (hae) bav -ar f 


Qui dA. 
3B = 


With regard to the total energy EP" introduced in (5.62), the last equation can be recast in 
the form 


ee | 


Tapi dA. (5.80) 
8,B 2 
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5.4. Numerical investigations 


Comparison with (5.61) confirms that the present scheme correctly reproduces the balance 
of energy. Note that the fulfillment of the directionality condition (5.69) is of paramount 
importance for obtaining an energy consistent scheme. 


Eventually, we aim at the discrete balance of entropy. To this end, multiply (5.73)3 by 
(D,n). and take into account (5.66)3, to obtain 


ab Tbn¢1 — Tbn _ h d 
(D-n)aH CAT. I MGE " VN*dV 
1 [E 
- | V | = |- Kł VOR dV — ? dA. 
B a) * S 0,8 Oligo 


The last equation can be rewritten as 


Lore) + Don" i (Chyi - CD) av 
B 


1 1 Jayi 
=at fy Ha Rh | av at | ? dA. 
B (s) l * E Oligo 0,8 Oligo 
eS 


>0 


The discrete derivatives of the entropy density, namely D-n" and Don", respectively, 
satisfy an associated directionality condition in analogy to (5.69). Accordingly, the left- 
hand side of the last equation equals SE — S. where S? is the total entropy of the 
discrete system introduced in (5.64). Comparison with (5.63) confirms that the present 
scheme correctly reproduces the balance of entropy. 


5.4. Numerical investigations 


As has been shown in the last section the newly developed schemes are thermodynam- 
ically consistent and comply with the balance laws for linear and angular momentum, 
respectively. These structure-preserving properties hold for arbitrarily large time-steps. 
Depending on the specific choice of the thermodynamic state variable r € {@,n, u}, 
we obtain three alternative EME schemes. The respective scheme is denoted as (EME),. 
For example, choosing the temperature as state variable leads to the temperature-based 
(EME)g scheme. 


We apply the newly developed (EME), schemes to representative numerical examples 
dealing with finite-strain thermo-elastodynamics. In particular, we treat the three examples 
which have been previously investigated in Section 4.5. In Chapter 4 the time discretization 
was confined to the mid-point rule which led to a genuine lack of numerical stability. In 
particular the lack of stability of the mid-point rule is especially pronounced in the case 
of large time steps. Therefore, in the numerical examples presented below we focus on 
large time steps to show that the new (EME); schemes remain stable. 
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5.4.1. Flying L-shaped block 


First we consider the example given in Section 4.5.1 and choose the largest time-step size 
for all three (EME), schemes. The total linear momentum of the block is a conserved 
quantity since in the initial loading phase the external forces are equilibrated. As after 
the loading phase (t > 5s) no external torque is acting on the block the total angular 
momentum is a conserved quantity as well. All (EME), integrators under consideration 
are capable to conserve both momentum maps (up to numerical round-off), independent 
of the time step size which can be observed from Fig. 5.4, where representative numerical 
results of the (EME), integrator are shown. 


40711 104 
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03 

5 

A 9 

= 

—2 — 

0 100 200 0 100 200 


— L; (EME),,-At = 0.4 —— J, (EME),,-At = 0.4 


—— L; (EME),,-At = 0.4 
— Ls (EME),,-At = 0.4 


— J, (EME),,-At = 0.4 
— J, (EME),,-At = 0.4 


Figure 5.4.: L-shaped block: Algorithmic conservation of linear momentum (EME),, scheme (left), Total discrete 
angular momentum (EME),, scheme (right) 


During the loading phase (0 < t < 5) the total energy of the block increases, whereas after 
the loading phase the total energy has to be a conserved quantity. All (EME), integrators 
exactly reproduce this conservation law (up to numerical round-off), see Fig. 5.5. This 
is in sharp contrast to the mid-point based schemes. In particular, as has been shown in 
Section 4.5.1, all mid-point based schemes exhibit numerical instability even though the 
formulation in terms of the internal energy density is capable to conserve the energy. 
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Figure 5.5.: L-shaped block: Total energy (EME), schemes (left), Incremental change of total energy (EME) 


schemes (right) 
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Figure 5.6.: L-shaped block: total entropy (EME), schemes (left), Incremental change of total entropy (EME)- 


schemes (right) 


Since the (EME), integrators are consistent with the second law of thermodynamics 
(independent of the time-step size), the total entropy is a non-decreasing function over 
time, as can be observed from Fig. 5.6. The above investigations indicate that, in contrast 
to the mid-point based schemes (see Section 4.5.1) the (EME), integrators are numerically 
stabile even for large time steps. After the loading phase, the environment of the present 
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example can be classified as thermally perfect in the sense of [146]. Therefore, £ introduced 
in (4.92) plays the role of a Lyapunov function. All of the (EME), integrators correctly 
deliver a decreasing value of £ over time (Fig. 5.7) thus confirming numerical stability. 
This is in sharp contrast to the mid-point based schemes dealt with in Section 4.5.1. 


104 
3 
E 
S 28 
2.6 ET 
0 20 40 60 80 100 120 140 160 180 200 220 240 
t [s] 
—— (EME),,-At = 0.4 
—— (EME),-At = 0.4 


—— (EME), -At = 0.4 


Figure 5.7.: L-shaped block: Lyapunov function computed with the (EME), schemes 


Eventually, the motion of the L-shaped block is illustrated in Fig. 5.8 with snapshots at 
successive points in time. In addition to that, the distribution of the temperature over the 
block is shown. 


360.00 


336.67 


313.33 


290.00 


Figure 5.8.: L-shaped block: Snapshots of the motion along with the temperature distribution over the block at 
t € {0, 32, 64, 96, 128, 160, 192, 224}, obtained with the (EME)g scheme and time step At = 0.4s 
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5.4. Numerical investigations 


5.4.2. Rotating disc 


Next we consider the example given in Section 4.5.2 and choose the largest time-step 
size for all three (EME), schemes. As neither external loads (0,8 = ()) nor displacement 
boundary conditions (0,,B = Ø) are acting on the disc, the thermomechanical system at 
hand has translational and rotational symmetry. Therefore the corresponding momentum 
maps are first integrals of the motion. All (EME), integrators at hand are capable to 
conserve the respective momentum map. Representative numerical results of the (EME);, 
integrator are shown in Fig. 5.9. During initial period t € [0, 4]s the total energy of the 
system is expected to increase due to the prescribed heat flow into the system. After 
initial period the system is classified as closed and therefore the total energy of the system 
should stay constant. All (EME), integrators are capable to correctly reproduce the first 
law of thermodynamics (see Fig. 5.10) in contrast to the mid-point based schemes, where 
only the formulation in the internal energy density was in accordance with the first law of 
thermodynamics, see Section 4.5.2. The total entropy of the system is expected to increase 
due to the prescribed heat flow into the system during the initial period t € [0, 4]s. As 
the system is closed after the initial period, the total entropy of the system should be a 
non-decreasing function, whereby the irreversibility is caused by heat conduction. 


.10713 .101 


— L, (EME), -At — 0.1 ——] (EME), -At =0.1 
— L; (EME),-At = 0.1 — J; (EME),-At = 0.1 


Figure 5.9.: Rotating disc: Total linear momentum (left) and total angular momentum (right) 
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— (EME), -At = 0.1 — (EME),-At = 0.1 


Figure 5.10.: Rotating disc: Total energy (left) and incremental change of total energy (right) 


All (EME), integrators correctly reproduce the second law of thermodynamics as can 
be observed from Fig. 5.11. Again this is in contrast to the mid-point based schemes 
investigated in Section 4.5.2, where only the formulation in terms of the entropy density 


was shown to be consistent with the second law. 
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Figure 5.11.: Rotating disc: Total entropy (left) and incremental change of total entropy (right) 
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5.4. Numerical investigations 


To investigate the numerical stability of the (EME), schemes further, we consider again 
the Lyapunov function £ introduced in (4.92). After the initial period the system is closed 


and therefore £ has to decrease with time. 
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Figure 5.12.: Rotating disc: Lyapunov function (EME), schemes 


All (EME), integrators correctly reproduce this behaviour as can be observed from Fig. 5.12. 
This again confirms numerical stability of the EME schemes at hand and is in contrast to 


the mid-point based schemes investigated in Section 4.5.2. 


Figure 5.13.: Rotating disc: Error in the position (left) and error in the velocity (right) 
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Figure 5.14.: Rotating disc: Error in the absolute temperature 


Further we investigate the order of accuracy of the present structure-preserving methods. 
Therefore we introduce the relative La error norm in the position, the velocity and the 
absolute temperature, respectively 


lec? . lvy [IO ul. 
*^ lel ' 7- lv Sl. ' 


FE av] 


In the above equations e, with e € {y, v, O} is the reference solution at time t = 4.02s 
calculated with the smallest time step size. We consider the motion of the rotating disk in 
the time interval [4, 4.02]s. To this end, the simulation is started at time t— 4s by using the 
data obtained with the (EME)g scheme and At— 0.04s. As can be observed from Figs. 5.13 
and 5.14 all schemes exhibit second order accuracy in y, v and ©. Eventually, the motion 
of the disc is illustrated in Fig. 5.15 with snapshots at successive points in time. In addition 
to that, the distribution of the temperature over the disc is shown. 


tol 


I| e llc; 
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360.00 


336.67 


313.33 


290.00 


Figure 5.15.: Rotating disc: Snapshots of the motion at successive points in time t € (0, 4, 8, 12, 16, 18, 24, 
28}s, and corresponding temperature distribution, calculated with the (EME)g scheme and At = 0.1s 


5.4.3. Rotating disc in a thermally perfect environment 


Finally we examine the example given in Section 4.5.3 and consider only the largest 
time-step size. Again we focus on the temperature-based (EME)e integrator, which makes 
possible to impose temperature Dirichlet boundary conditions in a standard manner. As in 
the previous example the system at hand has both translational and rotational symmetry. 
Accordingly, the corresponding momentum maps are exactly conserved by the (EME)g 
integrator, see Fig. 5.16. 
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Figure 5.16.: Rotating disc in a thermally perfect environment: Total linear momentum (left), total angular 
momentum (right) 
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According to the temperature boundary conditions heat is flowing out of the disc. The 
corresponding time-evolution of the total energy is depicted in Fig. 5.17 (left). In addition 
to that, the total entropy is as well decreasing over time, see Fig. 5.17 (right). 
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Figure 5.17.: Rotating disc in a thermally perfect environment: Total energy (EME)g scheme (left), and total 
entropy (right) 


By design of the thermally perfect environment, functional £ introduced in (4.92) is a 
Lyapunov function throughout the whole simulation time. Correspondingly, the (EME)o 
integrator yields a steadily decreasing value of £ even for very large time steps (see Fig. 
5.18). In particular, in contrast to the mid-point rule (cf. Section 4.5.3), the present scheme 
does stay numerically stable. 
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Figure 5.18.: Rotating disc in a thermally perfect environment: Lyapunov function 
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5.4. Numerical investigations 


Eventually, both the motion and the evolution of the temperature distribution are shown 
in Fig. 5.19 with snapshots at successive points in time. 
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Figure 5.19.: Rotating disc in a thermally perfect environment: Snapshots of the motion at t € (0, 10, 20, 30, 40, 
50, 60, 70) s and corresponding temperature distribution, calculated with the (EME)g scheme and At = 0.0875s 
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6. GENERIC-based numerical 
methods for large-strain 
thermo-viscoelasticity! 


Large-strain thermo-viscoelasticity is described in the framework of GENERIC. To this end, 
a new material representation of the inelastic part of the dissipative bracket is proposed. 
The bracket form of GENERIC generates the governing equations for large-strain thermo- 
viscoelasticity including the nonlinear evolution law for the internal variables associated 
with inelastic deformations. The GENERIC formalism facilitates the free choice of the 
thermodynamic variable. In particular, one may choose (i) the internal energy density, (ii) 
the entropy density, or (iii) the absolute temperature as the thermodynamic state variable. 
A mixed finite element method is proposed for the discretization in space which preserves 
the GENERIC form of the resulting semi-discrete evolution equations. The GENERIC- 
consistent space discretization makes the design of structure-preserving time-stepping 
schemes possible. The mid-point type discretization in time yields three alternative 
schemes. Depending on the specific choice of the thermodynamic variable, these schemes 
are shown to be partially structure-preserving. Numerical investigations are presented 
which confirm the theoretical findings and shed light on the numerical stability of the 
newly developed schemes. 


6.1. GENERIC-based formulation of large strain 
thermo-viscoelasticity 


The present Chapter relies on a material description of large-strain thermo-viscoelasticity 
within the framework of GENERIC. The developments presented herein extend the work 
in the previous Chapters 4 and 5 on thermo-elasticity to the realm of inelastic deformations. 
In the GENERIC framework the time-evolution of any functional A is governed by 


dA 


a TE} [^s]. (6.1) 


1 This Chapter is based on [4] 
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The evolution equation (6.1) is valid for closed systems in which the boundaries are 
disregarded. We consider a continuum body with material points X = X ^e4 in the 
reference configuration B C IR? (Fig. 6.1). The material coordinates X ^ refer to canonical 
base vectors e4 € R. Here and in the sequel, the summation convention applies to 
repeated indices. 


ej 


&3 
Figure 6.1.: Reference configuration B and deformed configuration p(B, t) at time t. For now the focus is on 


the motion of isolated thermomechanically coupled solids. That is, external tractions acting on the boundary as 
well heat fluxes across the boundary are disregarded. 


Within the Lagrangian description of continuum mechanics the deformed configuration 
of the body at time t is characterized by the deformation map q : B x T +> R3, where 
I = [0,T] is the time interval of interest. Accordingly, the placement at time t of 
the material point X € B is given by x = q(X,t). The corresponding velocity field 
v: B x Z R? is defined by v = q, where a superposed dot denotes the material time 
derivative. The deformation gradient is given by 


F = dx, (6.2) 


or, in components, (F)', = Ox'/OX4. Furthermore, the right Cauchy-Green tensor 
reads 


C=F'F. (6.3) 
In the GENERIC (6.1), we consider functionals of the form 
A=A(p,B;T, C;') = no C;)aV, (6.4) 
B 
where p = pv is the linear momentum density and p : B +> R+ is the mass density in 


the reference configuration. Moreover, 7 : B x Z ++ Ris a generalized thermodynamic 
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6.1. GENERIC-based formulation of large strain thermo-viscoelasticity 


field which may be chosen from among three alternatives, 7 € {6,7,u}. The three 
alternative fields are the absolute temperature 6, the entropy density 7, and the internal 
energy density u. To describe inelastic deformations, we make use of the internal variable 


C-1: Bx T. R3X3, Since we focus on isotropic inelastic deformations, C7 ! is assumed 
p P p 


to be symmetric, i.e. C = C; . Therefore the density function a in (6.4) is subject to 
the same isotropic restrictions as the Helmholtz free energy function, to be introduced in 
Section 6.2. 


GENERIC (6.1) decouples reversible and irreversible processes. In particular, the Poisson 
bracket b j is responsible for reversible processes, while the dissipative (or friction) 
bracket |, -] embodies irreversible processes. Accordingly, these two brackets constitute 
fundamental entities of the GENERIC framework and shall be specified next. Since inelastic 
deformations are purely irreversible, the Poisson bracket retains its thermo-elastic form 
(see eq. (4.34)) 


{A,B} -/ (5,0 + Div (san) ) - pb — (5,0 + Div (Fr) - ópa dV , 


(6.5) 


where A and B are arbitrary functionals of the form (6.4). Furthermore, the functional 
derivatives in (6.5) are given by 


ôpa = Oya — DivOga , 


dpa = Opa, 
d n (6.6) 
ô-a = 0,a, 
ðc-14 = c-a . 
In (6.5), 7 denotes the entropy density giving rise to the total entropy of the system 
S(p, T, C7) = n GC, dV. (6.7) 
B 


The dissipative bracket for the thermo-viscoelastic problem at hand can be additively 
decomposed into a part due to heat conduction and a part taking into account inelastic 
deformations: 

[A.B] = [A.B]... + 48] (6.8) 


cond inel ` 


In analogy to thermo-elasticity with heat conduction, the dissipative bracket accounting 
for heat conduction is given by (see eq. (4.35)) 


u 67a 2 6,5 
[A, B] ona = fy (=) -0°KV (55) dV. (6.9) 


133 


6. GENERIC-based numerical methods for large-strain thermo-viscoelasticity 


Here, K = K(C,r) is the material conductivity tensor, for which we assume that 
K” = Kanda: K : a > 0 holds for all a € R?. An important element of the 
GENERIC framework is that the absolute temperature takes the form 


= O,u(F, T, Cj!) 
O(F, T, C, 1) = nE, 6. . (6.10) 
25 25? p 


In (6.9) and (6.10), u denotes the internal energy density which, together with the kinetic 
energy and the potential of dead loads, constitutes the total energy of the system 


HE 1. s 
Epro - | & !p:puF.r C!) - be) dV. (6.11) 


Here, b : B ++ R? represent prescribed body forces which, for simplicity, are assumed to 
be dead loads. 


Concerning the contribution of inelastic deformations, we introduce the following form 
of the dissipative bracket (see Section 6.2 for further details): 


[A.B] 


inel 


ôa k 2 örb a = 
- E (Zooi; z 2) ON :2 xe 50-100, Jj av. 


(6.12) 


The fourth-order inelastic material flow tensor N has the properties NT =N (major 
symmetry) and M : N : M > 0 (positive semi-definiteness) for all M given in (6.17). 


It can be easily verified that the Poisson bracket (6.5) is skew-symmetric ((.A4, B) = 
—{B, A}), the dissipative bracket (6.12) is symmetric ([.A, B] = [B, A]) and further sat- 
isfies the non-negativity condition [A, A] > 0. Moreover, the two brackets satisfy the 
fundamental degeneracy (or non-interaction) conditions {A,S} = 0 and [A, E] = 0. In 
conjunction with these properties of the two brackets, GENERIC (6.1) ensures for closed 
systems (i) the conservation of total energy (d£ /dt = 0), and (ii) a non-decreasing total 
entropy (dS /dt > 0). Due to these structural properties, the GENERIC-based formulation 
offers an ideal starting point for the design of structure-preserving numerical methods. 


6.1.1. Local evolution equations 


We provide a short outline of the local evolution equations emanating from GENERIC (6.1). 
With regard to (6.1) the total energy € generates the reversible contribution through the 
Poisson bracket {A, €}. Using expression (6.11) for the total energy along with formulas 
(6.6) for the variational derivatives, Poisson bracket (6.5) yields 


{A, E} = | (a+ Div (Zen) -p ‘p+(b+DivP)-dpadV, — (6.13) 
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where the first Piola-Kirchhoff stress tensor P — P(F, 7, 6; takes the form 
P = öru — OFN . (6.14) 


The irreversible contribution to GENERIC (6.1) is generated by the total entropy S through 
the dissipative bracket [A, S]. Inserting expression (6.7) for the total entropy into the 
dissipative bracket (6.8) leads to 


[A,S] = fy (=) Qa + [2 (52e 6; -30,100;") :N:MdV, 
B T B T s 5 
(6.15) 


where the material heat flux vector Q = Q(F, 7, c,') is given by 
2 1 
Q=0 KV = -KVO, (6.16) 


and M = M(F, 7, C 1) denotes the material representation of the Mandel stress tensor, 
which takes the form 


M = 2 (5. iu — 6dg-10) C7 (6.17) 


It is important to realize that in (6.14), (6.16) and (6.17), O is the temperature field which has 
been introduced in (6.10). These relationships are representative for the GENERIC-based 
formulation. 


With regard to the left-hand side of GENERIC (6.1) and expression (6.4) for functional .A, 
the time derivative of A can be written as 


d LL 
al, (50: + a Bt brat dea: c dv. (6.18) 
Substituting (6.18), (6.13) and (6.15) into GENERIC (6.1), we arrive at 
0 = f 5e: (e — p p)dV 
B 


«f dpa: (p — (b+ DivP)) dV 
B 


5, (6.19) 


. 2 -1Y. : _p; 970 nl 
+ [5.0 (+- F (opc; JVM) Div (Sen) p pdV 


5, aa a 
n (zz) Qaa: (cz «2w Mpeg?) dV. 
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This equation has to hold for arbitrary functionals A of the form (6.4). Standard arguments 
now lead to the local evolution equations emanating from GENERIC (6.1): 


p=p'p, 

p=b+DivP, 

. 2 AY Lug. 1 zd Los (6.20) 
T= Ju (86:6; ) N :M ga : V(p ! p) 3, rA. 


C; --2(N:M)CG;!. 


We stress again that in the above formulas the first Piola-Kirchhoff stress tensor P, the 
Mandel stress tensor M and the heat flux vector Q are given by formulas (6.14), (6.17) 
and (6.16), respectively. 


6.1.2. GENERIC-based weak form of the IBVP 


Starting from the GENERIC-based local evolution equations (6.20), we deduce the weak 
form of the initial boundary value problem (IBVP) pertaining to finite-strain thermo- 
viscoelasticity. To this end, we decompose the boundary OB of the continuum body 
into a displacement boundary 0,,8, on which ọ = Ẹ, and a traction boundary 0,B, on 
which PN = t, where @ and t are prescribed functions for t > 0 (Fig. 6.2). Moreover, 
0,B U 0,B = OB and 0,8 0,B = 0. 


ƏB 


ej 
ea 


Figure 6.2.: Mechanical part of the IBVP. Note that t — PIN denotes prescribed external Piola tractions acting 
on the current boundary expressed per unit area of the reference boundary ôs B. 


Similarly, for the thermal part we consider the subsets of the boundary ô, and 0,8, with 
the properties 0-6 U 0,B = OB and 0,B 0,B = () (Fig. 6.3). The thermodynamic state 
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variable is prescribed on 0,B, i.e. T = 7, whereas the heat flux is prescribed on Ö,B, i.e. 
Q- N =. Both 7 and q are assumed to be given for t > 0. 


p 
N A N 


I 


ea 


Figure 6.3.: Thermal part of the IBVP. Note that g = Q - N is the prescribed rate of heat transfer across the 
current boundary expressed per unit area of the reference boundary 0,B. 


We aim at the determination of the motion e : Bx Z > R3, the linear momentum 
p: B xT R?, the thermodynamic state variable 7 : B x Z > R, and the inelastic 
deformation Cj, 1: B x e R?*?, The unknown fields are subject to initial conditions 
of the form y(-,0) = X, p(-,0) = pVo, T(-,0) = ri?! and C, (., 0) = I in B. Here, Vo 
is a prescribed material velocity field, and 7"! is a prescribed field of the thermodynamic 
state variable 7 € {0, n, u}. The unknown fields are determined by applying a space-time 
discretization to the weak form of the IBVP at hand. 


The weak form of the IBVP can be obtained by scalar multiplying the local evolution 
equations (6.20) by suitable test functions and subsequently integrating over domain B. 
The standard procedure yields 


0= [wo (o - v) dV 


w w 
-V = ).QdV "qdA 
|, (55) = + 4B Anu" 


+ | wep : (G57 «20 iM) e!) dV , 
B p 
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where the velocity field v = p~'p has been introduced. In weak form (6.21), Wy, Wp : 
Bo RÌ, w, : B > Rand Wo: Be R?*° are test functions that have to satisfy the 
p 


boundary conditions w = 0 and wp = 0 on 0, B, and w, = 0 on 9, B. 


6.1.3. Balance laws 


We verify the pertinent balance laws of the IBVP at hand. For that purpose it suffices 
to consider the pure Neumann problem (ie. 0,B = 0,B = OB). We start with the 
total linear momentum of the continuum body defined by L = f g P dV. Choosing 
(Wy, Wp, Wr, Wez) = (0, £,0,0), where £ € R? is arbitrary but constant, weak form 


(6.21) leads to 
E = 
— = bdV tdA). 6.22 
erae). ea 


Due to the arbitrariness of € € IR?, (6.22) coincides with the balance law for linear 
momentum. Note that the parentheses on the right-hand side of (6.22) contain the resultant 
external force exerted on the continuum body (see also Fig. 6.2). 


The total angular momentum relative to the origin of the inertial frame is defined by 
J= Ip? x p dV. Choosing (Wy, Wp, Wr, Wa-ı) = (p X €,€ x », 0,0), weak form 


(6.21) yields 
dJ E 
E (fexpav+ [ oxa). (6.23) 


Note that the symmetry condition FP" = PF” has been employed. Since € € R? is 
arbitrary, (6.23) corresponds to the balance of angular momentum. The parentheses on 
the right-hand side of (6.23) contain the resultant external torque about the origin (see 
also Fig. 6.2). 


Next, we substitute (Wo, Wp, Wr, Wo, 1) = (-b,p^ 1p, Ou, dc; ıu) into weak form 


(6.21). A straightforward calculation leads to the balance law for total energy 


d 1 s 
so 'p-p+u av= [b-oav+ f (o'p-t-q) dA. (6.24) 
dt B aB 


For the balance of entropy we choose (Wy, Wp, Wr, Wc-:) = (0,9 'p, 0-7, Og-1) in 
p p 
the weak form (6.21) to obtain 


-f e : F + ðn? + 20¢-10 : (N :M) C,!) —2 (86:7) iN: M) dV 


+ [ 2, am: 6; ep eti qav+ f Onn dA 
ru op Oru 


a ES ET : M NU 
4G 5M:N:M- v(3) e'xv (5) ) ave f. uaa. 
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Here, use has been made of formula (6.10) for the temperature along with expressions 
(6.16) for the material heat flux vector and (6.17) for the Mandel stress tensor, respectively. 
The above equation can be rewritten as 


di —gdA- | V(—)-O?KV(—)+=M:N:MdV>0, (625 
dt * Jas 6° iE (s) (5)*s N 20, (625) 
en 
=Deond= 0 =Dinei 20 


where Dine! is the local production of entropy due to inelastic deformations and Dona is 
the local production of entropy due to heat conduction. The last equation corresponds to 
the Clausius-Duhem form of the second law of thermodynamics (see, for example, [136, 
Sec. 5]). 


6.2. Inelastic part of the dissipative bracket 


The present model for large-strain thermo-viscoelasticity relies on the introduction of 
the internal variable C, 1: Bx T o R3*3, whose time-evolution accounts for local 


inelastic deformations. Following [148], C,, ! can be associated with the multiplicative 
decomposition of the deformation gradient [149] 


F=F-.F,, (6.26) 


into an elastic part F. and an inelastic (or viscous) part Fp. Decomposition (6.26) gives 
rise to the relationship 


Cp = FF). (6.27) 
The restriction to the isotropic case implies that the free energy takes the separable form 
(see, for example, [150]) 
V = V(F,7,C,') = V*(F,7) + V"*(F, r, C), (6.28) 
where W°4 is the equilibrium part and V"*3 is the non-equilibrium part. 


Since inelastic deformations are purely irreversible in nature, they lead to a contribution to 
the dissipative bracket in GENERIC (6.1), cf. (6.8). To derive the inelastic dissipative bracket 
(6.12), we resort to the dissipative bracket derived in [99] within a temperature-based 
framework for GENERIC. In particular, [99] consider functionals of the form 


K(o,p,,F,) = / 8(o, F, p,6, F,) dV , (6.29) 
B 


with associated density function a(¢, F, p, 0, Fp). The inelastic dissipative bracket from 
[99] is given by 


[A,B]. =| (o a J (ON: (aos usd B) dV, (6.30) 
ine. B p P Oot p P 
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where u(F, 0, F,,) is the internal energy density and N is a fourth-order inelastic flow 
tensor which has the properties NT = N (major symmetry) and A : N : A > 0 for all 
A € R?*?, In components, these properties read 


(NI LENS. dado QUAIN DOSE 0s (6.31) 


As has been shown in [99], inelastic dissipative bracket (6.30) comes along with the flow 
rule 


F,=-N:>, 


where 


X= Of, — 00r, : 


In the last equation, 7((F, 0, F,,) is the entropy density of the temperature-based formula- 
tion. 


6.2.1. Change of variables 


We perform a change of variables to transform the inelastic bracket (6.30) to the present 
setting which is based on functionals of the form (6.4). In particular, to link the current 
density functions a(q, F, p, 7, cz") to those in (6.29), we express the generalized ther- 
modynamic state variable 7 € (0,1, u} in terms of the temperature by inverting relation 
(6.10) to get 


r=7(F,0,05'): (6.32) 
Moreover, relation (6.27) implies 
-1 _ p-lp-T 
C, =F, Fp 3 (6.33) 
Now, the two density functions under consideration can be connected through 


ap, F, P; 0, F,) == a(p, F, p,7(F, 0, Cy), C,') , 


where relationships (6.32) and (6.33) are employed on the right-hand side of the last 
equation. A straightforward application of the chain rule to the last equation yields 


a = O,aO0gT ; 


i B (6.34) 
Op, à = Op, a + O,a0p,T . 


Note that 69a = 09a, and df, à = Op, à. Furthermore, with regard to (6.10) and (6.32) we 


have 0 = O(F,7(F, 0, c Co) from which follows that 


1 = 090 = 8,9097 , 
0= Og, 9 = 0, OOF, 7 + Og, 0 ! 
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We thus obtain 


SENI 

idm ch (6.35) 

Jr- Og, O i 
P 8,0 
Substituting from (6.35) into (6.34) yields 
o-a 
Opa = SE : 
a &.a (6.36) 
Op, à = OF, a — PX oF © : 


Note that 6-a = Ö,a and óp, a = Op, a. Making use of (6.36), the terms in the parenthesis 
of (6.30) can be rewritten as 


Sn ee Op a= LN u—Óóp,a. (6.37) 
Oo u 


Taking into account the relationship 
=F -1 
r,a = —2F, 0o :aC, ; (6.38) 
the inelastic dissipative bracket (6.30) can be recast in the form 


[A.B], 


inel 


òa —1 = 6,b =1 ESI 
= E (Fee; ues = ôcz10C, ) : ON : 2 (22e; 6; — 66,100, dV. 


This bracket coincides with the one introduced in (6.12). Note that in the above formula, 
óc -14 = Og-1a. Moreover, N is the material form of the fourth-order flow tensor N in 
p 


(6.30). In components, 
(N) Ts (6.39) 


Accordingly, the material flow tensor inherits symmetry and positive semi-definiteness. 
That is, in components, 


(N);5;2(N)É;; and (M) (NYL (M)g" 20 (6.40) 


for all M introduced in eq. (6.17). As has been shown in Section 6.1.1, the above inelastic 
dissipative bracket gives rise to flow rules of the form (cf. (6.20)4) 


6, =-2(N:M)C,'. (6.41) 
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It is worth noting that this flow rule can be viewed as material version of the viscoelastic 
evolution equation 


Lobe = 20 vp be , (6.42) 


derived in [148, 150] and being applied in, e.g., [151]. In the last equation, be = F.F; = 


FC; FT, Lybe = FC; FT, qued 2FOüc WFT, and Y^! is an isotropic, positive 
definite fourth-order tensor. Using these relationships, flow rule (6.42) can be recast in 
the form 


C; =-2F-! (v^ : (F-"Mr7)) FC;!, (6.43) 


where the relation 


2COcW'* = 205.9" 30,1 =M, 


has been employed. Note that in the last equation, definition (6.17) of the Mandel stress 
tensor has been taken into account. Comparing (6.43) with present flow rule (6.41), leads 
to the conclusion that the respective fourth-order flow tensors are related by 


(NJE ed yap ae qu bue. 


We further remark that the present viscoelastic evolution equation (6.41) can also be 
brought into the form 
C, -2€0,N : M. 


This version of the viscoelastic evolution equation has been used in [106], see also [76]. 
Thus, we conclude that the newly proposed inelastic dissipative bracket (6.12) gives rise 
to evolution equations for the internal variables that have been previously developed in 
the context of finite deformation thermo-viscoelasticity. 


6.3. Discretization in space 


Concerning the space discretization of the present GENERIC-based formulation, we essen- 
tially apply an isoparametric finite element approach (see, for example, [145]). Our main 
goal is to achieve a GENERIC-consistent space discretization in the sense of Section 2.1.2. 
In particular, à GENERIC-consistent space discretization ensures that the discrete formu- 
lation inherits the fundamental balance laws for both the energy and the entropy from 
the underlying continuous formulation (cf. Section 6.1.3). 
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We first restate the governing equations of the IBVP to be discretized in space and time. 
With regard to the GENERIC-based weak form (6.21), we consider the following set of 
equations: 


(6.44) 


2 1 
0= fw. (4-2 (2¢,uc;") Ns M+ “de Vv) dV 
B Nr 


-fv (=) av + f WradAdV. 
B Ur 3B Ur 


Here, the first equation represents the local form of the kinematic relationship % = v, 
while the second equation is the viscoelastic evolution equation. The third and fourth 
equation serve the purpose to introduce the new fields u, and 1j. This procedure facilitates 
a mixed finite element approach which turns out to be crucial for a GENERIC-consistent 
discrete formulation. Due to the arbitrariness of the test functions wy, and w;, , (6.44)3 
and (6.44), impose the conditions u, = 0,u and n, = +n, respectively. Moreover, 


P= Opu — OFN, 


E (6.45) 
M =2 (ðs; 1u- 68.0) C; 

are the GENERIC-specific representations of the first Piola-Kirchhoff stress tensor and the 
material Mandel stress tensor previously introduced in (6.14) and (6.17), respectively. Sim- 
ilarly, the material heat flux vector Q has been introduced in (6.16). We emphasize again 
that the GENERIC-based formulation is based on expression (6.10) for the temperature 
field. In the present mixed formulation this implies 


gol. (6.46) 
N 


The finite element method is based on finite dimensional approximations of the following 
quantities 


Q^ = N*q, , 
vh = Nv, (6.47) 
7 = Nz, , 
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and 
h a 
uz = N"(ur)a ; 
4 b r)a (6.48) 
7, = N*(n-)a 
As before, the summation convention applies, where a = 1,..., N, and N denotes the 


total number of nodes in the finite element mesh. Moreover, N*? : B — Rare the nodal 
shape functions with associated nodal values q,, v, € R? and T4, (ur)a, (Nr)a € R. 
Analogous approximations are used for the test functions Wu, , Wn, , Wp, and w;, denoted 


h ch yh h 
by Wa, Wn, Wy, and w7. 


In what follows, we summarize the space-discrete version of (6.44). Nodal collocation of 
kinematic equation (6.44); yields 


da — Va; (6.49) 


for a = 1,...,N. Viscoelastic evolution equation (6.44)2 is collocated at the integra- 
tion points X, € B used for the numerical evaluation of the volume integrals in (6.44). 
Accordingly, 


(6,5, = -2(N,:M3)(C;),, (6.50) 


for g — 1,..., G, where G denotes the total number of integration points. Here and in 
the sequel, index g indicates evaluation at the integration point X, € B. For example, 


(C; ), = Qr (Xat) ’ 


7 (6.51) 
My = 2 (8: - 950: n;) (C5 s 


where 
ug=u (Fa (Cs) ; 


(6.52) 
Ng = n (Fi. oe (o) ’ 


are the internal energy and entropy densities at point X,. Furthermore, the discrete 


deformation gradient at X4, F7, and the discrete generalized thermal variable 75. follow 
from (6.47),,3 and thus take the form 


rh = q, (1) ES VN*(X ) , 
: 1 (6.53) 


iP N*(X4)r«(t) . 


Similarly, the discrete temperature at X, follows from (6.46) and is given by 


9, - (ut), (6.54) 
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where, in view of interpolations (6.48), (u®), = ub(X,,t) and (gh), = nb(X,,t). In 
the discrete setting, the fields u? and rj? are determined through (6.44)3 4. In particular, 
inserting the approximations (6.48) along with the corresponding formulas for wh ; wr 
into (6.44)3 4, we obtain 


g=1 
K (6.55) 
0= 5 N*(X,) (cn, — N'(X,)(u.);) Wg , 
g=1 
for a = 1,...,N. To calculate the spatial integrals, appropriate numerical quadrature 


formulas of the form 
G 
f f(X)dV & M f(Xj)u, , (6.56) 
B g=1 


have been applied to obtain (6.55)?. Here, wg play the role of generalized weighting 
coefficients resulting from the specific quadrature rule along with the isoparametric 
description of reference domain B. Now, (6.55); can be solved for the nodal quantities 
(u,)4, a = 1,..., N. To this end, we introduce the components H“ of the positive definite 
Gram matrix [H*'], 


H* = V N*(X,)N*(X,)u, , (6.57) 
g=1 
so that (6.55); can be rewritten as 
G 
H*^(u,)y = V ^ N*(X,)0-ugu, . (6.58) 
g=1 
The components Ha, of the inverse Gram matrix, [H.»] = [H^^|-!, satisfy the relation- 
ship 
Ha H” = 6;° , (6.59) 


where ô ° denotes the Kronecker delta. Now, interpolation formula (6.48); along with 
(6.58) lead to the result 


ub = N*Hay >, N’ (X), ugw, à (6.60) 


? The summation over g will always be stated explicitly, so that the summation convention does not apply to 
index g. 
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which will be utilized below. Similarly, interpolation formula (6.48)2 in conjunction with 
(6.55)2 yield the result 


G 
ne = N° Has X N*(Xg)ðrNgWwg ; (6.61) 


g=1 


Next, we turn to the discretization of (6.44)5, which can be done in a straightforward way 
to obtain 


G G 
M'^v, — S  N*(X,)bw, + 3 P, VN"(X;)u, 20. (6.62) 
g—1 g—1 


For simplicity, we have neglected the contribution of the external tractions which could 
be easily added to the above equation. In the above equation, first Piola-Kirchhoff stress 
tensor P, at point X, € B is given by 


P, = Fug — OgOF Ny . (6.63) 
Moreover, in (6.62), the components M^ of the mass matrix are given by 
G 
M? = V5 o(X,)N*(X,)N"(X,)u, . (6.64) 
g=1 


We further introduce nodal momentum vectors p^ conjugate to nodal position vectors qa 
through the standard relation 


p^ = M**y, . (6.65) 


Eventually, we consider the space-discrete version of (6.44). Straight-forward application 
of our approach yields 


S V2N*(X 
0 =H = p e (86: u4 (€; a) N, : Mw; 
ed | H N (6.66) 
"E _ 9^ 8gn.V N*(CX. EP V (2) 2 : 
Ye Du E DV |“ Gaus) Qus 


where the material heat flux vector Q, at point X, € B is given by 


sene (t) 


For simplicity, in the space-discrete evolution equation (6.66) for the generalized nodal 
thermal variable 7), the term accounting for heat transfer across the boundary has been 
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neglected. To summarize, the resulting evolution equations for the space-discrete system 
at hand can be written in the form 


Ga =Ma p” ’ 
G G 
p^ =) N*(Xj)b,u, - M P,VN"(X,)u, , 
g=1 g—1 


G 
: 2N°(X,) E 
Ta =Hab ) Bc (85: us(C5 Jo) : Ng : Mywg 
g=1 v 


G G 
Y N*(X : Y N*(X 

— Ve: Hab PU a) opn, YN (X5)w, + Aap V ( = 2) Q, wg , 
gl T79 T/g 


(6,5, 2 2 (N, : M5) (C53), , 


(6.67) 


for 1,..., N and 1,..., G. In (6.67), May stands for the components of the inverse mass 
matrix satisfying Map M"* = 6;*. The set of equations (6.67) constitute nonlinear first- 
order ordinary differential equations for the determination of the unknowns which can 
be collected in the state vector 


1 N -1,4B —1\4B 
$e (di: Titae TN (Cp Jj? yag e jo ) (6.68) 
In the sequel, state vector (6.68) will be viewed as column vector. In particular, this implies 
that the six independent components (0708, g = 1,...,G, of the internal variable 
(C, 2 g (at quadrature point X.) are arranged in a column vector. 


The set of evolution equations (6.67) fits into the GENERIC framework for discrete systems, 
as shown next. 


6.3.1. GENERIC-consistent space discretization 


Our discretization approach presented above is GENERIC-consistent in the sense of [2] 
and thus can be framed in the context of GENERIC for discrete systems. To see this, the 
set of evolution equations (6.67) needs to be put into the form 


z= LV£(z) + MVS(z). (6.69) 


Here, the focus is again on closed systems in which the boundary contributions are 
disregarded. In analogy to (6.1), the time-evolution of the state vector is decomposed 
additively into a reversible part generated by the total energy € and an irreversible part 
generated by the total entropy S. Poisson matrix L needs be skew-symmetric, while 
friction matrix M. needs be symmetric positive semi-definite. 
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In the above equation, the total energy of the discrete system under consideration is given 


by 
G 
ER i i 
£(z) = 5 Map“ p° + Y^ [u (Fo, 73, (CRo) - a N*QG)bOXS)| wy - (6:70) 


This is the space-discrete version of total energy (6.11). Similarly, the space-discrete 
version of total entropy (6.7) takes the form 


S) - 37» (F; g? faal Ce) Wg - (6.71) 


g—1 


To get the gradient of the total energy, V£(z), we consider the derivative of (6.70) with 
respect to time. Accordingly, the left-hand side of (6.70) yields 


G 
d AB e 
ape (2) = 04,6 da Ops E H° Or, + WED NG rave (©, er (6.72) 
while the right-hand side of (6.70) gives 
G 
d -a b . a a 
m = p“ - P Ma cd: 5 [drug V. N* (Xg) — N*(X;)b(X,)] wg 
= (6.73) 
esae g)Ó ugwg + 1585 c; iin ACE us 


g—1 


Comparing (6.72) with (6.73) yields the following expressions for the respective derivatives 
of the total energy 


G 

Bau = Y ru V N* (X) — N'(X;)b(X,)] w9 , 
g=1 

Op E = VL 


(6.74) 
0,,E = Lra, J)O rUgWyg , 


In particular, inserting from (6.74)3 into (6.60), we obtain the important relationship 


= N*HaO,£. (6.75) 
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Similarly, the derivatives of total entropy (6.71) take the form 
G 


04,8 = 5 Og, VN" (Xu); , 
g—1 
Op. S =0, 
G (6.76) 
0,,5 = SON (XK )0.ngw, , 
g—1 
Os + Ong 
aj, AB g AB 
a(r), 9(C;) 


Inserting from (6.76)3 into (6.61), we obtain 
= N*°H»0,,8 . (6.77) 
Now, guided by discrete GENERIC (6.69), the evolution equations in (6.67) can be recast. 
In particular, taking into account (6.74)2, kinematic relation (6.67); can be rewritten as 
à, = Op. (6.78) 


Next, evolution equation (6.67)2 for the nodal momentum vectors together with (6.63) and 
(6.54) yields 


G G h 
=- Y eu TNX) - N*OG)bQG wy Y EEI VN" Kye 
g=1 g 


Taking into account (6.74) and (6.75), we obtain 


G 


N*(X 
p° = —O4,E + OnE Hoc X zh a) oen, VN" (Xj) (6.79) 
T79 


Concerning the evolution of the nodal thermal variable Ta, (6.67)3 along with (6.51)2, (6.54) 
and (6.77) lead to 


G b c 
: ANXN X) 
fa Hoo) Tab) 


g=1 


(86: (C75) : Ws : (8e (075) wy Ha 
= DX M odds TE (86: u,(€7),) Wg 
— OpeE - Ha 


NOX N*(X 
+ Ha 089 ( ED aeo C 2) toy Habs 


ge (ub), (u7 )g 


2 oeny NX g)Wg 


(6.80) 
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Eventually, evolution equation (6.67)4 for the internal variable, together with (6.51)o, (6.54) 
and (6.77) result in 


(6,5, E NM (N, : (852: u5(054)) (C, Hadr,S (6.81) 


t 40, e : (86: n,(C7),)) (€, ), . 


Next, we aim at introducing relation (6.76)4 into the last term of (6.81). To this end, we 
introduce the fourth-order tensor M with components 


(MNAR = (C Te BAD EOS" 3 (6.82) 


Note that M enjoys major symmetry, (M)4?!/ = (M)!7AB, due to the major symme- 
try of the material flow tensor N and the symmetry of C,.. Now, (6.81) can be recast in 
index form 


"a Ou N?(X o Os 
(C; d - —4(M) 381 = B C o) 07,8 4 4— (A)2PT7 Adj 
(C5!) (ms Wg 8(C;!), 


(6.83) 


Altogether, evolution equations (6.67) pertaining to the state variables of the discrete 
system at hand can be recast in the form 


Ga ES pE ’ 
pP" =— ôa E +1, On€ , 


G 
3 ] Os 
Tag = (aye Op» E + Mab 07,5 + so aCe , (6.84) 
g=1 P jg 
(C, 77 = mj OS c mz?! ZR 
Cp’), 
where 
G c 
1°, e? hy Ogr V Now, , 
g—1 TEY 
G c d 
N N 
Mab = Hac > Q2V (ae) -K,V | —*+ | wg Ha 
2.9 (cay, ) e | Gay, ) v 
G c 
+H Y ANON? Qus (M)ABII Jug lips (6.85) 
ac h h _1, AB g ayes 49 ? 
g=1 (uana. O(Cp 2 O(Cp 1 
Ou N2 
AB ABIJ g g 
mAB = —A(MDAPH C9 9 er, 
DR * aces)! Qs 
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6.3. Discretization in space 


For simplicity of exposition, in (6.85), N = N“(X,,). Note that the properties Mab = Mba 


and A = more hold. 


Now, evolution equations (6.84) give rise to specific forms of Poisson matrix L and friction 
matrix M in GENERIC (6.69). In particular, the Poisson matrix takes the form 


0 I 0 0 
I 0 [*;] 0 

L= ; (6.86) 
0 pu 0 0 
0 0 0 0 


where I is the identity matrix (with appropriate dimension corresponding to the partition- 
ing of the state vector (6.68)), and matrix [17;] consists of vectors 1°, defined in (6.85). 
Specifically, we have 


D eed, 


Bela: GR a (6.87) 


It can be observed that Poisson matrix (6.86) is skew-symmetric. Furthermore, the friction 
matrix is given by 


0 0 0 0 

Se 0 0 0 0 l (6.88) 
0 0 [ma] | [m4] 
0 0 [ma2]" | [m45] 


Here, the block matrices [mas] € RN*N, [m7] € RN*96, and [m} #17] contain the 
components defined in (6.85)2—4 and take the form 


Mi c0 MIN [mii] n [mg] 
mal=| : oo. 0: fs ads : : 
MN1 CS MNN [mi] wee [m& x] 
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and 
[maBH] ... 0 
0 es [más] 


Note that the last matrix is block-diagonal and symmetric. It can be easily observed that 
friction matrix (6.88) is symmetric and positive-semidefinite. 


6.3.2. Conservation properties 


We next verify the conservation properties of the semi-discrete formulation of the closed 
system dealt with in Section 6.3.1. Since the evolution equations pertaining to the semi- 
discrete formulation can be brought into GENERIC form (6.69), (i) conservation of total 
energy, and (ii) non-decreasing total entropy are automatically satisfied. To see this, we 
first consider the time-derivative of the total energy, 


Tel) =VE(z)z 
= VE(z)LVE(z) + VE(z)MVS(z) , 


(6.89) 


where (6.69) has been used. The first term on the right-hand side of the last equation 
vanishes due to the skew-symmetry of Poisson matrix (6.86). The second term vanishes 
too, since the non-interaction condition 


MV£(z) =0. (6.90) 


holds. The last equation can be easily verified by a straight-forward calculation. Accord- 
ingly, (6.89) yields the conservation law d£ /dt = 0. Similarly, the time-derivative of the 
total entropy yields 


(6.91) 
= VS(z)LV£(z) + VS(z)MVS(z). 


It can be verified by a straight-forward calculation that the second non-interaction condi- 
tion 


LVS(z) =0. (6.92) 


is satisfied. In addition to that, the positive semi-definiteness of friction matrix (6.88) 
implies the result dS /dt > 0. In particular, this result represents the semi-discrete version 
of (6.25) (apart from the boundary term in (6.25) originating from heat flux across the 
boundary). Thus, the total entropy of the closed system at hand is non-decreasing, due to 
the irreversible nature of heat conduction and visco-elastic deformations. 


152 


6.3. Discretization in space 


Since the material response is assumed to be frame-indifferent (or objective), specific 
symmetry properties are inherent to the discrete system under consideration. In particular, 
invariance under rigid rotations implies satisfaction of the following conditions (see 
Appendix C.1 for further details): 


0=q, X ôa E, 
da is (6.93) 
0=q, x E ’ 
for all b = 1,..., N. We tacitly assume that no resultant external torque is acting on the 


system (i.e. the right-hand side of (6.23) is assumed to vanish). The semi-discrete version 
of the angular momentum relative to the origin is given by 


j= / pe^ x v dv = qa x p^, (6.94) 
B 


where formulas (6.47)1,» along with definition (6.65) of the nodal momentum vectors p° 
have been used. The time-derivative of (6.94) reads 


d 
`J” =å, xp*4 x p“ 
dt ER SP (6.95) 


= Map x p* qa x( dt IE) , 


where (6.67), and (6.84)? have been used. Now the right-hand side of the last equation 
vanishes due to (i) the symmetry of Ma» together with the skew-symmetry of the vector 
cross product, and (ii) symmetry conditions (6.93). Result dJ" /dt = 0 implies conservation 
of total angular momentum. 


6.3.3. Choice of the thermodynamic state variable 
We shortly outline the impact of the specific choice of the thermodynamic state variable, 
T € (0,9, u}, on the structure of GENERIC (6.69). For simplicity, in this section we neglect 


body forces, that is, b = 0 and still focus on closed systems. 


Choosing the internal energy density as thermodynamic state variable, ie. 7 = u, the 
total energy (6.70) takes a particularly simple form given by 


G 
5 1 a b h 
E(z) = 5 Map": p taste. (6.96) 
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As before, ub stands for u^ (X,, t). In particular, interpolation formula (6.53)2 gives rise 


to ub = N*(X,)u,(t). The derivatives of the total energy in (6.74) simplify to 


04,8 —0, 
3p- Ê = Map” ’ 
R G 
8,6 = M N*(Xj)u,, (6.97) 
g—1 
aÊ 
L0. 
-1,AB 
(C; 3s 
Consequently, 
d'a , 
—€(z) = VE(z)2 


G (6.98) 
= Map" : p^ + 3 5 N*(X,)üsu, . 
g=1 


Moreover, for the choice 7 = u friction matrix (6.88) attains a particularly simple block- 
diagonal form, since map = 0, and coefficients m.» only contain contributions due to 


heat conduction (cf. (6.85)). 


Choosing the internal entropy density as thermodynamic state variable, i.e. 7 = 1j, yields 
a particularly simple form of the total entropy (6.71) given by 


G 
S(z) =X nw, , (6.99) 
g—1 


where interpolation formula (6.53)2 gives rise to n} = N“(X,)n.(t). Consequently, the 
derivatives of the total entropy in (6.76) simplify to 


04,8 = 0, 
IpaS =0, 
, G 
On, = M N*(X,)u, , (6.100) 
g—1 
as 
=). 
—14AB 
9(C; jo 
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Thus 
di _ 
a (z) 2 VS(z)z 
G (6.101) 
= 5 N? (X )nawg 
g—1 


Moreover, the choice 7 = 77 leads to 1”, = 0 (see (6.85)), so that Poisson matrix (6.86) 
yields a particularly simple form. 


In contrast to the above considerations, selecting the total temperature as thermodynamic 
state variable, i.e. 7 = 0, essentially does not lead to any simplifications. We eventually 
remark that these conclusions also affect the discretization in time, which will be treated 
next. 


6.4. Discretization in time 


We now turn to the discretization in time of the semi-discrete GENERIC-consistent evo- 
lution equations derived in Section 6.3.1. To this end, we focus on a representative time 
interval [tn, £4.41] with corresponding time-step size At = t441 — tn. We aim at second- 
order accurate, implicit time-stepping schemes based on the mid-point rule. Application 
of the mid-point rule to (6.69) yields 


Zn+1 — Zn = Atb(z,41)V£(2541) + AtM(z, 1) V5S(z,,1) - (6.102) 
Here, Z» stands for the discrete vector of state variables at time tn, and 


Zu+3 = (Zn + Zn41) - 


NI = 


Provided that the state variables z,, are given, the state variables z,,;1 can be determined 
by solving (6.102). 


6.4.1. Partially structure-preserving schemes 
Next, we check whether, or under what conditions, structure-preserving properties hold 
in the discrete setting. In this connection, we shall see that the specific choice of the 


thermodynamic state variable r € (0,1, u} plays a crucial role. It can be easily verified 
that non-interaction conditions (6.90) and (6.92) are still satisfied in the sense that 


0 
6.103 
6 (6.103) 
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To see whether the fundamental balance laws are correctly reproduced in the discrete 
formulation, we proceed along the lines of the time-continuous formulation in Section 
6.3.2. In particular, concerning the balance of energy, similar to (6.89), we consider 


V£(2541) (Zaà1 — Zn) = AIVE (2,41) i L(z,, 1) V£(2,41) 
+ AtVE(z,41) -M(z141)VS(zu41) ; 


[vr 


where (6.102) has been used. In analogy to the time-continuous case the right-hand side of 
the above equation vanishes due to the skew-symmetry of L(z, , i ) and non-interaction 
condition (6.103);. Thus 


V£(2,41) (541 — Zn) 70. 
On the other side, 
V£(254.1) ` (Zu+1 m Zn) * E(Zn41) = E(Zn) (6.104) 


in general. This inequality complies with the well-known fact that the mid-point rule 
is not capable to conserve nonlinear first integrals in general. However, there exists the 
exceptional case related to the choice 7 = u, for which (6.104) turns into an equality. 
This is due to the fact that for 7 = u the total energy takes the form (6.96) and thus É(z) 
is merely quadratic. Since the mid-point rule preserves quadratic first integrals (see [6, 
Sec. 4.4.2]), the choice 7 = u yields a structure-preserving scheme which is capable to 
conserve total energy. 


Concerning the evolution of total entropy, guided by (6.91), we consider 
VS(z,,1) x (Zn+1 = Zn) m AtVS(Z,41) i L(z,41)VE(2n42) 
+ AtVS (2,41) : M(z,,1)VS(2,,1) ] 


where again (6.102) has been used. Employing non-interaction condition (6.103)2, we 
obtain 


VS(z,,1) - (Zn+1 — Zn) = AtVS (2,41) : M(2,42)VS(2n42) 20, (6.105) 


where the positive semi-definiteness of friction matrix M(z, ,. i ) has been taken as a basis. 
On the other hand, in analogy to (6.104), we have 


VS(z,41) ‘(Zn41 — Zn) Z S(Zn41) — S(Zn) - (6.106) 


This implies that, despite the encouraging result (6.105), the mid-point scheme in general 
does not guarantee a non-decreasing entropy. However, there again is an exception related 
to the choice 7 = n. For this particular case, the total entropy takes the form (6.99) and 
thus S(z) is merely linear. Accordingly, the choice 7 = n turns inequality (6.106) into 
an equality and the entropy-based mid-point scheme is therefore capable to correctly 
reproduce the second law of thermodynamics in the discrete setting. 
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We eventually verify that all mid-point-based schemes under consideration are capable to 
conserve angular momentum. The incremental change of angular momentum (6.94) can 
be written in the form 


h h 
Jasi Jn = da; X (Pati = Phi) — Paya X (dios - da, ) : (6.107) 
Mid-point scheme (6.102) gives rise to 
b 
Qa, = da, = AtMavP y+ ’ 


a (6.108) 
Pnai — Papi = At (-9u.& 2143) + 1% (Zn42) Om&(2u+4)) : 


Inserting from (6.108) into (6.107) yields 


Jh E J = Atda, 1 x (-d4.€ (#44) + 1% (Zn44) OnE Zn44)) _ AtMavPi | x p, 


nri * 
(6.109) 
Symmetry conditions (6.93) imply 
0= da. 1 x Oa, € (Zn+4) ’ 
DS he (6.110) 
0= a x I") (2n44) 5 
2 


Inserting from (6.110) into (6.109) and taking into account the symmetry of Ma» together 
with the skew-symmetry of the vector product leads to the result ds = Jh. 


To summarize, depending on the choice for the thermodynamic state variable we get three 
alternative mid-point schemes which are partially structure-preserving. Correspondingly, 
the resulting schemes are denoted by (EM),, (energy-momentum scheme related to 7 = u), 
(ME),, (momentum-entropy scheme associated with 7 = 7), and Mg (momentum scheme 
related to 7 — 0). 


6.5. Numerical investigations 


In this section, the alternative mid-point type schemes newly developed in the present 
Chapter are applied to representative numerical examples dealing with finite-strain thermo- 
viscoelastodynamics. Depending on the specific choice for the thermodynamic state 
variable r € (0,9, u}, the following methods are applied: 


u | (EM), | Energy-Momentum consistent scheme 
n | (ME), | Momentum-Entropy consistent scheme 
0 | (M)a Momentum consistent scheme 


A|! 
Il 
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Following the previous Chapters we shall focus on momentum maps associated with 
symmetries of the mechanical system at hand, and the balance laws associated with 
the two fundamental laws of thermodynamics. In this connection we also consider the 
functional 


L=E-09S$. (6.111) 
According to [146], for certain types of environments, £ is a natural Lyapunov function and 
thus qualifies as estimate for the numerical stability of the schemes under consideration. 


In each time step, the schemes emanating from (6.102) generate a system of nonlinear 
algebraic equations for the determination of the state variables z,+1. To this end, we 
apply a Multilevel-Newton algorithm, see [152] and references therein for more details. 


6.5.1. Material model 


In order to particularize the Helmholtz free energy density (6.28) used in the numerical 
examples we start from a temperature-based description. In particular, we consider 


—neq 


P(C, 0, C5!) =P (C, 0) +0" (C,9,Cz!), 
Y (C, 0) = 41(C) + v»(8) — (0 — %)v3(J) , (6.112) 
(C, 0, C;") == V1 vise (C; C;") , 


— neq 


vj 


where 
E i 2 2 
w(C) =3 C:1-3-2lgJ - 2 (J-— 1) t Wya(J) ; 


e - 2 
V4 visc(C, 0,7) = ^ (c : C, j en 3 i 2 log Je m 3 (Je va n?) T Wyolvise( Je) , 


V3(8) = c (0 — 6 — 61og(0/89)) , 
vs(J) = 3BWyaa(J) , 


à+? 
Wal) = — (log J)? + (J —1)2) , 
re + Site 
Wvolwise(Je) = n ((log Je)? + (Je = 1) ) S 


(6.113) 


Here, J = vdetC is the determinant of the deformation gradient and Je = 
4/ det (CC; ). Further u, He, A and X. are prescribed parameters, c > 0 is the spe- 


cific heat capacity at constant deformation, f is the coefficient of thermal expansion, and 
0o is the reference temperature. We refer to [1] and the references therein for a detailed 
investigation ofthe thermoelastic part of the specific Helmholtz free energy density (6.112). 
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Further, for the viscoelastic part we refer to [74] and the references therein. For simplicity 
we assume incompressible material behavior. Quasi incompressible material models for 
finite-strain thermo-viscoelasticy are considered in, e.g., [142, 153]. It is now a straightfor- 
ward exercise to calculate further quantities emanating from (6.112), depending on the 
specific choice for the thermodynamic variable 7 € {u, 0, n} (see also [1] for additional 
details). In particular, the temperature-based formulation yields 


7(C, 0, C; !) = clog(8/60) + vs(J) , 
u(C, 0, C,') = V1(C) + Wrvise(C, C;") + c(8 v 0o) + bows (J) E 


The formulation based on the entropy density gives 


WC, n, C, ) =n, 
ee E EE n—w3(J) 
&(C,n, C5!) = vn (C) + Pryise(C, C; EC ae 60) TENOR 


Moreover, the formulation based on the internal energy density leads to 


— (C) — Yı ,vise(C, CZ") sa) + vs(J) 


U 
n(C, u, C, !) = clog ( + 
cho C 
&(C,u, C; !) =u. 


Concerning the constitutive equation for the material heat flux vector, we assume thermally 
isotropic material, with material conductivity tensor given by 


K-kJC |. (6.114) 


Here, k is a prescribed coefficient of thermal conductivity and, as before, J = ‚/det(C). 
Finally, the constitutive equation for the isotropic fourth-order material inelastic flow 
tensor is given by (see [69] or [148] and references therein for the spatial representation 
of the isotropic fourth-order inelastic flow tensor) 


1 1 1 
- ——[[C^!oI| C- <IeI ——IG«I 6.115 
(| ol 181) * ger, (6.115) 


where vp > Qand vy > 0 represent the deviatoric and volumetric viscosities, respectively, 


and where (A © BJ? 57? = 3 [(A)?€ (B) 5? + (A)22(B) 5°]. 


6.5.2. Flying L-shaped block 


The first numerical example deals with the L-shaped block depicted in Fig. 6.4. 
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Figure 6.4.: L-shaped block: Discretised block with initial temperature distribution and mechanical boundary 
conditions (left), load function over time (right) 


The spatial discretization of the block relies on 117 tri-linear finite elements leading to 
224 nodes. The initial temperature field is varying linearly over the height (x3 direction) 
of the block. In particular, at x3 = 0, the initial temperature is prescribed as 0,, while at 
x3 = h, the temperature is prescribed as 0,. The whole block is assumed to be thermally 
insulated (q = 0 on 0,8 = OD). Starting at rest, Piola traction vectors ta and t, are acting 
on two parts of the boundary surface of the block (Fig. 6.4). The external loads are applied 
in the form of a hat function over time. In particular, the traction vectors are given by 


256/9 t for Os<t< 2s, 
ta = —t, = f(t) | 512/9 | Pa, f@)=<4-t for 2s<t<4s, (6.116) 
768/9 0 for t>4s. 


Table 6.1 provides a summary of the data used in the simulations. During the loading 
phase (t < 4s) the time step size for all simulations is At = 0.05, such that all systems 
start from the same energy-and entropy level directly after vanishing external loads. No 
Dirichlet boundary conditions are applied. Since in the initial loading phase the external 
forces are equilibrated, the total linear momentum of the block is a conserved quantity. 
In addition to that, after the loading phase (t > 4s) no external torque is acting on the 
block. 
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Lamé parameters 


Specific heat capacity 
Expansion coefficient 
Thermal conductivity 
Viscosities 


Ref. temperature 
Mass density 
Initial temperature 


Geometry 


Newton tolerance globoal 
Newton tolerance local 
Simulation duration 
Time step 


At 


997.5 
4544 
49.875 
272.2 
100 
2.233 - 1074 
10 

500 
100 
293.15 
100 


1078 

107° 

300 

0.05, 0.5, 0.6 


Pa 

Pa 

Pa 

Pa 
JK^!m-3 
Kk 

WK !m^! 
Jsm~? 
Jsm~? 

K 

kgm ? 

K 

K 

m 

m 


S 
S 


Geometry 


Table 6.1.: L-shaped block: Data used in the simulations 


Consequently, the total angular momentum is conserved as well. All of the integrators un- 


der consideration are capable to exactly conserve both momentum maps (up to numerical 
round-off), independent of the chosen time step. This can be observed from Fig. 6.5, where 
representative numerical results of the (EM),, integrator are shown. After the loading 
phase the total energy must be a conserved quantity. As expected, the (EM),, scheme 


does exactly reproduce this conservation law (up to numerical round-off), see Fig. 6.6. In 
contrary, the schemes (M)g and (ME), are not capable of conserving the total energy for 
larger time step sizes. Depending on the time step size, both schemes tend to increase the 


energy which can be observed from Fig. 6.7. Typically, such energy blow-ups eventually 
lead to a failure of the iterative (Newton-Raphson) solution procedure. In the diagrams 
the break down of the simulation is indicated by vertical lines. 
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404 


J; [Nsm] 
oO 


— L, (EM), -At = 0.05 —— J, (EM), -At = 0.05 
——L, (EM) -At = 0.05 — J, (EM), -At = 0.05 
— Ls (EM) -At = 0.05 — J; (EM),,-At = 0.05 


Figure 6.5.: L-shaped block: Algorithmic conservation of linear momentum (EM), scheme (left), Total discrete 
angular momentum (EM), scheme (right) 


0 100 200 300 


— (EM),,-At = 0.05 — (EM),-At = 0.05 


— (EM),,-At = 0.5 


— (EM),-At = 0.5 
— (EM), -At = 0.6 


— (EM), -At = 0.6 


Figure 6.6.: L-shaped block: Total energy (EM), scheme (left), Incremental change of total energy (EM). scheme 
(right) 
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-105 -105 


—— (M)g-At = 0.05 — (ME),- At = 0.05 
— (M)g-At = 0.5 — (ME),-At = 0.5 
— (M) -At = 0.6 — (ME),-At = 0.6 


Figure 6.7.: L-shaped block: Total energy (M)g scheme (left), Total energy (ME), scheme (right) 


Regardless of the capability of the (EM). scheme to conserve the total energy, the simu- 
lation still breaks down at about the same point in time as the break down of (M)g and 
(ME), occurs. The numerical instability of the (EM), is accompanied by a nonphysical 
decrease of the total entropy as can be observed from Fig. 6.8. In fact, the total entropy 
ought to be a non-decreasing function of time. In contrast to (EM), and (M)e, (ME), is 
capable to correctly adhere to the second law of thermodynamics, independent of the time 


step (Fig. 6.9). 
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— (M)g-At = 0.05 


— (EM),-At = 0.05 
— (M),-At = 0.5 


— (EM), -At = 0.5 


— (EM), -At = 0.6 — (M),-At = 0.6 


Figure 6.8.: L-shaped block: Total entropy (EM),, scheme (left), Total entropy (M)g scheme (right) 
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Indeed, in each time step the total entropy does increase, as can be observed from Fig. 6.9. 
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— (ME), -At = 0.05 
— (ME),-At = 0.5 
— (ME),-At = 0.6 


— (ME),-At = 0.05 
— (ME),-At = 0.5 


— (ME),-At = 0.6 


Figure 6.9.: L-shaped block: total entropy (ME); scheme (left), Incremental change of total entropy (ME); 
scheme (right) 


For sufficient small time step sizes, the incremental change of total entropy is governed 
by 
Sad = [hai Pos. „av + al Die „av (6.117) 


Bere mer due to Deona Contribution due to Dine) 


for all mid-point based schemes. The two contributions to the discrete evolution of the 
incremental change of total entropy are visualized in Fig. 6.10 and Fig. 6.11. 
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—— (EM),-At = 0.05 
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Figure 6.10.: L-shaped block: Contributions to incremental change of entropy (EM). scheme (left), Contributions 
to incremental change of entropy (M)g scheme (right) 
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Figure 6.11.: L-shaped block: Contributions to incremental change of entropy (ME); scheme 


After the loading phase, the environment of the present example can be characterized 
as thermally perfect in the sense of [146]. Thus £ defined in (6.111) plays the role of 
a Lyapunov function that has to decrease with time. However, the partially structure- 
preserving schemes (EM),,, (ME), and (M)g do not correctly reproduce this behavior, as 
can be seen from Figs. 6.12 and 6.13. That is, depending on the time step and the duration 
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of the simulation, all of the schemes inevitably exhibit numerical instabilities characterized 


by increasing values of £. 


— (M)g-At = 0.05 
— (M)e-At = 0.5 
— (M)g-At = 0.6 


—— (EM),,-At = 0.05 
— (EM),-At = 0.5 
—— (EM), -At = 0.6 


Figure 6.12.: L-shaped block: Lyapunov function (EM),, scheme (left), Lyapunov function (M)g scheme (right) 


— (ME), -At = 0.05 
— (ME),-At — 0.5 
— (ME),-At = 0.6 


Figure 6.13.: L-shaped block: Lyapunov function (ME); scheme 


Eventually, the motion of the L-shaped block is illustrated in Fig. 6.14 with snapshots at 
successive points in time. In addition to that, the distribution of the temperature over the 


block is shown. 
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360.00 


333.33 


306.67 


280.00 


Figure 6.14.: L-shaped block: Snapshots of the motion along with the temperature distribution over the block at 
t € (0, 40, 80, 120, 160, 200, 240, 280] s, obtained with the (EM), scheme and time step At = 0.05s 


6.5.3. Rotating disc 


The second example deals with a rotating disc subjected to prescribed heat flow over part 
of the boundary surface (Fig. 6.15). The spatial discretization of the disc is based on 200 
tri-linear finite elements leading to a total of 360 nodes. 


wo 1 
0o 
ea 0.5 
R” 
ei 0 
02 4 6 810 
t [s] 
q 


— ft) 


Figure 6.15.: Rotating disc: Initial configuration and thermal boundary conditions (left), function f (t) for the 
prescribed heat flow over part of the boundary surface (right) 


The initial velocity distribution over the disc results from a prescribed angular velocity 
wo € R? and is given by 


AE 
vo(X) = Wo X X, Wo = 1 Sond 
1 S 
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The initial temperature of the disc is homogeneously distributed and equal to the reference 
temperature 0p. In an initial period of time, t € [0, 4]s, heat flow is prescribed over one 
quarter of the lateral boundary surface (Fig. 6.15). In particular, the heat flow into the disc 
is described by 


2000W 


71m? 


sin(Zt) for 0<t<4s, 


0 for t>4s. 


q=- 


A plot of function f(t) can be found in Fig. 6.15. The rest of the boundary surface of the 
disc is assumed to be thermally insulated (q = 0). Note that the prescribed heat flow 
vanishes after t = 4s. For t > 4s, the complete boundary surface of the disc is assumed 
to be thermally insulated (q = 0 on 9;B = OB). Then the environment of the disc is 
thermally perfect in the sense of [146]. A summary of the data used in the simulations 
of the rotating disc can be found in Table 6.2. During the loading phase (t < 4s) the 
time step size for all simulations is At = 0.04s, such that all systems start from the same 
energy-and entropy level directly after vanishing external loads. 


Material parameters À 3000 Pa Geometry 
H 750 Pa 
He 120 Pa 
A 37.5 Pa 
Specific heat capacity c 150 JK !m-? 
Expansion coefficient B ivi x i PR st 
Thermal conductivity k 20 WK m^! N J} 
Viscosities Vp 50 Jsm-^? 
Vy 10 Jsm-? 
Ref. temperature 0o 293.15 K 
Mass density p 893 kgm ? 
Radius r 08 m 
R 2 m 
Thickness t 04 m 
Newton tolerance global £g 1075 - 
Newton tolerance local El 107? - 
Simulation time T 30 s 
Time step At 0.04, 0.08,0.1 s 


Table 6.2.: Rotating disc: Data used in the simulations 


Due to the fact that neither external loads act on the disc, nor displacement boundary 
conditions are imposed, the mechanical system at hand has translational and rotational 
symmetry. Consequently, the corresponding momentum maps are first integrals of the mo- 
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tion. All integrators under consideration are capable to conserve the respective momentum 
map. Representative numerical results are shown in Fig. 6.16. 


— L, (ME),-At = 0.04 ——Jı (MB),-At = 0.04 
— L, (ME),,-At = 0.04 — J, (MB),-At = 0.04 
— L; (ME),,-At = 0.04 — J; (MB),-At = 0.04 


Figure 6.16.: Rotating disc: Algorithmic conservation of linear momentum (EM),, scheme (left), Total discrete 
angular momentum (EM); scheme (right) 


Since heat flow into the system is prescribed in the initial time period [0, 4]s, the total 
energy is expected to increase. For t > 4s, the system is closed and the total energy should 
stay constant. Again the (EM),, scheme is capable to correctly reproduce the first law 
(Fig. 6.17). 
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—— (EM),,-At = 0.04 
—— (EM),,-At = 0.08 
— (EM),-At = 0.1 


—— (EM),,-At = 0.04 
— (EM),,-At = 0.08 


— (EM),-At = 0.1 


Figure 6.17.: Rotating disc: Total energy (EM). scheme (left), Incremental change of total energy (EM). scheme 
(right) 


However, despite the algorithmic energy conservation (for t > 4s), the (EM),, scheme is 
not devoid of numerical instabilities, depending on the time step. The corresponding point 
in time of the break down of the simulation is indicated with a vertical line in the diagrams. 
At about the same points in time, (ME),, and (M)e break down as well (Fig. 6.18). For these 
schemes the break down is accompanied by a sudden increase of the total energy leading 
to the divergence of the Newton-Raphson iterations. For the considered duration of the 
simulation (t € [0, 30]s), a time step of At = 0.04s is small enough to retain numerical 
stability of the three partially structure-preserving schemes at hand. 
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10? 10? 


—— (M)g-At = 0.04 — (ME),-At = 0.04 
— (M)g-At = 0.08 — (ME),-At = 0.08 
—— (M)g-At = 0.1 — (ME),-At — 0.1 


Figure 6.18.: Rotating disc: Total energy (M)g scheme (left), Total energy (ME); scheme (right) 


Due to the prescribed heat flow into the disc, the total entropy of the disc is expected to 
increase in the initial time period [0, 4]s. For t > 4s, the system is closed and the total 
entropy should be a non-decreasing function of time. The (EM),, scheme does not correctly 
reproduce the second law as can be seen from Fig. 6.19. Accordingly, the divergence of 
the iterative solution procedure is accompanied by a nonphysical decrease of the total 
entropy. The (M)e closely adheres to the second law as can be observed from Fig. 6.19. 
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— (M)g-At = 0.04 
— (M)g-At = 0.08 
— (M),-At = 0.1 


—— (EM),,-At = 0.04 
—— (EM),,-At = 0.08 
— (EM),-At = 0.1 


Figure 6.19.: Rotating disc: Total entropy (EM), scheme (left), Total entropy (M)g scheme (right) 
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As expected, the (ME),, scheme is capable to exactly satisfy the second law, independent of 
the time step. This can be observed from Fig. 6.20. In particular, Fig. 6.20 (right) confirms 
that the change per time step of the total entropy is always positive. 
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— (ME),- At = 0.04 
— (ME),-At = 0.08 


— (ME),- At = 0.04 
— (ME),-At — 0.08 


— (ME),-At — 0.1 — (ME),-At = 0.1 


Figure 6.20.: Rotating disc: total entropy (ME); scheme (left), Incremental change of total entropy (ME); 
scheme (right) 


In addition, the two contributions to the discrete evolution of the incremental change 
of total entropy are visualized in Fig. 6.21 and Fig. 6.22 for the mid-point based schemes. 
Most of the contribution is due to conduction of heat, only about 5.9% is due to inelastic 


deformations in the given example. 
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mu Contribution due to Di; mum Contribution due to Dine 


Figure 6.21.: Rotating disc: Contributions to incremental change of entropy (EM). scheme (left), Contributions 
to incremental change of entropy (M)g scheme (right) 
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Figure 6.22.: Rotating disc: Contributions to incremental change of entropy (ME), scheme 


To shed further light on the numerical stability of the present schemes, we consider the 
Lyapunov function £ defined in (6.111). For t > 4s the system is closed and the function 
£ should decrease with time. As expected, the partially structure-preserving schemes 
(EM),,, (ME), and (M)g can not guarantee to correctly reproduce this behavior, depending 
on the size of the time step and the duration of the simulation (Figs. 6.23 and 6.24). In 
particular, it can be seen that the smallest time step, At — 0.04s, yields a stable numerical 
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simulation, at least in the considered time interval [0, 30]s. However, for larger time steps 
At = 0.08s and At = 0.1s, the numerical instability of each scheme becomes visible 


through the increase of £. 
10? 


.102 
2 4 
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— (M)g-At = 0.04 
— (M)g-At = 0.08 


—— (EM), -At = 0.04 
—— (EM), -At = 0.08 
— (EM),-At = 0.1 


Figure 6.23.: Rotating disc: Lyapunov function (EM), scheme (left), Lyapunov function (M)g scheme 


— (ME),-At = 0.04 
— (ME), -At = 0.08 
— (ME),-At = 0.1 


Figure 6.24.: Rotating disc: Lyapunov function (ME), scheme 


Eventually, the motion of the disc is illustrated in Fig. 6.25 with snapshots at successive 
points in time. In addition to that, the distribution of the temperature over the disc is 


shown. 
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Figure 6.25.: Rotating disc: Snapshots of the motion along with the temperature distribution over the block at 
t € {0, 4, 8, 12, 16, 20, 24, 28}s, obtained with the (ME), scheme and time step At = 0.04s 
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7.1. Summary! 


This thesis has addressed the formulation, analysis, implementation, and verification of 
second-order accurate structure-preserving schemes by utilizing the GENERIC structure. 


First, we have analyzed in Chapter 3 energy and momentum conserving schemes in 
the context of molecular dynamics. Conserving schemes have been employed for more 
than four decades for the accurate and robust approximation of Hamiltonian problems 
in mechanics. However, in the field of molecular dynamics, where time integrators are 
routinely employed and are key to their usefulness, the use of conserving schemes has 
barely been explored. These second-order implicit methods are an interesting alternative 
to other commonly used integrators, and we have proven that, in addition to exhibiting 
exact energy conservation, they are more robust than the midpoint rule, the canonical 
second-order implicit method. We focused on the design of energy-momentum (EM) 
schemes for molecular dynamics in the view of three issues that are characteristic and 
unique to these problems: (i) the numerical treatment of dynamics in periodic domains, 
(ii) the discretization of three-body potentials, and (iii) the study of interatomic functional 
potentials. None of these three topics had been previously studied in the context of 
conserving schemes, to the authors' knowledge. However, the three of them are key for 
their implementation and clearly have an impact on their performance. Some of the most 
practical results of this Chapter are new expressions for EM approximations in fairly 
general problems in molecular dynamics. These approximations account for the three 
key issues mentioned before and can be shown to preserve linear momentum and energy 
exactly, while exhibiting a remarkable robustness. 


In addition, in Chapter 4, we developed a GENERIC-based variational formulation of the 
IBVP for large-strain thermo-elastodynamics. The weak form of the problem directly 
emanates from the GENERIC evolution equation for open systems. This formulation 
makes possible the free choice of the thermodynamic state variable among the three 
options: (i) the absolute temperature 0, (ii) the entropy density 7, and (iii) the internal 
energy density u. The GENERIC-based weak form is particularly well suited for the design 
of structure-preserving numerical schemes. This observation already holds true for the 
discretization in time employing of the standard mid-point rule. Depending on the specific 


! The summary is based on [1-4] 
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choice of the thermodynamic state variable, specific balance laws are correctly reproduced 
in the discrete setting, independent of the size of the time-step. In particular, by using 
the internal energy density as the thermodynamic state variable, the resulting scheme 
is capable of consistently reproducing the balance of energy and thus has been termed 
energy-momentum (or (EM).) scheme. Similarly, choosing the entropy density as state 
variable leads to an entropy consistent method, termed momentum-entropy (or (ME),,) 
scheme. However, none of the three newly developed mid-point-based schemes does 
preserve all of the key balance laws under discretization. In particular, it was shown that 
numerical instabilities can occur if the time-step size is too large for a prescribed duration 
of the simulation. This leads to the conclusion that all major balance laws should be 
maintained in the discrete setting in order to enhance the numerical stability. 


Consequently, novel fully structure-preserving numerical methods for finite-strain ther- 
moelasticity with heat conduction were developed in Chapter 5, starting from the newly 
developed GENERIC-based weak form. The proposed EME schemes enable the free choice 
of the thermodynamic state in the same manner. Each choice of the thermodynamic 
variable 7 € {7,0,u} leads to a corresponding (EME), scheme. We have shown that 
two modifications of the mid-point-based approach are crucial for obtaining numerically 
stable EME schemes: (i) the introduction of specific projections in the wake of the space 
discretization and, (ii) the replacement of the standard derivatives of the internal energy 
density and the entropy density, respectively, by appropriate discrete derivatives. The in 
Chapter 2 introduced notion GENERIC-consistent space discretization is crucial with regard 
to the discretization in space. This notion asserts that a consistent space discretization 
should fit into the GENERIC framework for discrete systems. In particular, a GENERIC 
consistent space discretization automatically inherits the balance laws for energy and 
entropy, respectively, from the underlying continuous formulation. We have seen that 
the discretization in space, relying on standard Lagrangian shape functions, necessitates 
specific projections to reach a GENERIC-consistent space discretization. Even though the 
structure-preserving properties of the EME schemes are independent of the specific choice 
for T € (1,0, u}, the absolute temperature can be regarded as the preferred choice for the 
constitutive description. In addition to that, temperature Dirichlet boundary conditions 
can be applied in the standard manner by using the temperature-based formulation. It 
was shown that the EME schemes lead to a significant improvement in the numerical 
stability when compared to mid-point type schemes. 


Finally, in Chapter 6 the GENERIC-based approach was extended to more involved coupled 
thermomechanical problems, which also account for viscous (inelastic) deformations. 
Starting from a GENERIC-based formulation of large-strain thermo-viscoelasticity, we 
have developed alternative structure-preserving schemes based on the implicit mid-point 
rule. Depending on the choice of the thermodynamic variable, the plain mid-point rule 
yields partially structure-preserving schemes just as in the case for thermoelasticity. Also, 
in analogy to this case as well, these schemes cannot prevent numerical instabilities, as 
was shown in the numerical examples. 
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Upon summarizing the findings of the newly proposed GENERIC-based schemes, the 
following features are eminent: 


« The underlying GENERIC formulation leads to characteristic expressions such as 
those for the first and second Piola-Kirchhoff stress tensors, the Mandel stress 
tensor, and the absolute temperature. 


* The transformation properties of the GENERIC description facilitate the use of 
alternative thermodynamic variables such as the absolute temperature, the internal 
energy density, and the entropy density used in the present work. 


* The newly proposed material form of the inelastic dissipative bracket makes it 
possible to include viscoelastic effects into the GENERIC formulation, which are 
often used nonlinear evolution laws for the internal variables associated with 
inelastic deformations. 


e A GENERIC-consistent space discretization is essential when transitioning from 
infinitesimal dimensional to finite dimensional systems. This means that the evo- 
lution equations for the state variables of the semi-discrete system fit into the 
GENERIC framework for discrete systems. This is an essential prerequisite for the 
development of fully structure-preserving EME schemes. 


7.2. Outlook 


Based on the work contained within this thesis, some of the following lines of research 
are currently in progress or seem worthy of investing: 


« The GENERIC-based semi-discrete weak form for the thermo-viscoelastic problem 
is well suited for the design of structure-preserving schemes because the GENERIC 
structure is preserved after the discretization in space. Thus, starting from the 
GENERIC-consistent space discretization it is verly likely that EME numerical 
schemes with enhanced numerical stability and robustness will emerge when 
applying the discrete gradient operator in analogy to Chapters 4 and 5. 


e Volumetric constraints of incompressible rubber-like materials (e.g., polycarbonate) 
should be accounted for within the GENERIC framework, because viscous dissipa- 
tive features for such materials are usually associated with isochoric deformations, 
whereas the thermal coupling is naturally associated with volumetric deformations. 
The resulting GENERIC-based EME methods would greatly extend the scope of 
structure-preserving methods to thermodynamical systems with constraints. 


* The extension to non-smooth dissipative systems seems very promising, as effects 
like plasticity or damage are clearly irreversible and, therefore, the numerical 
method itself should align with this property. Incorporating an EME consistent 
time integration scheme would guarantee that the dissipated heat during such 
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irreversible transformations is not transferred to the mechanical field? A first 
attempt in this direction can be found in [124] for a thermo-elastoplastic spring, but 
the transcription to infinitesimal dimensional systems is, to the author's knowledge, 
yet to be made. 


Despite their good properties, structure-preserving numerical methods have not 
yet been implemented in common software packages; one reason for this is the very 
elaborate implementation of the classical discrete gradient operator. A promising 
attempt to circumvent this problem is the usage of an alternative discrete gradient 
operator introduced in [154] (see [155] for several contributions in this direction). 
This alternative version utilizes a polyconvex structure for the Helmholtz free 
energy density function which, paired with the tensor-cross product introduced 
in [156], leads to a very simple and elegant structure for the discrete gradient 
operator and could consequently makes structure-preserving numerical methods 
more appealing for commercial software packages. 


Interestingly, the class of " GENERIC-integrators" has recently been coined [157] and 
can be viewed as an extension of symplectic integrators for Hamiltonian systems to 
the dissipative regime. They may link the present work to the class of variational 
integrators. 


Finally, another interesting avenue for future study is the development of mixed 
finite elements within the GENERIC framework. The combination of GENERIC- 
based EME methods with mixed finite elements gives rise to even more stability 
and accuracy, as mixed finite elements are known for their improved performance 
compared to classical finite elements like the ones used in this thesis. Finally, a 
further improvement of the spatial resolution could be obtained by incorporating 
the framework of isogeometrical analysis [158]. 


? Solely energy-consistent schemes correctly balance the amount of energy contained within the system but are 
generally not concerned about the direction of energy transfer and thus do not generally satisfy the second 
law of thermodyamics (e.g., see results of (EM), methods in Sections 4.5 and 6.5). 
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A.1. Linear thermoelasticity 


We verify that the Helmholtz free energy density (4.36) yields an expression for the 
linearized stress tensor which is consistent with the theory of linear thermoelasticity. 
According to Section 4.2.4, the temperature-based expression for the first Piola-Kirchhoff 
stress tensor reads 

P = ygu — 00v 
Taking into account the functions 7(V, 0), and u(V qo, 0) summarized in Section 4.2.3, 
the first Piola-Kirchhoff stress tensor can be recast in the form 


P(F, 0) = Dui (F) — (0 — 60)5(J)cof(F) (A.1) 


It can be easily verified that the reference configuration is stress-free in the sense that 
P(L 05) — 0. To perform the linearization of the stress tensor (A.1) about the stress-free 
reference configuration, consider one-parameter families of the deformation gradient and 
the temperature of the form 

Fe =I+eVu 

0 = 0o c e(0 — 09) 


Here, the displacement field is denoted by u. Note that P(F°, 0?) = 0. The linearized 
stress tensor can now be calculated via 


d — 
o := —P(F^,0*) 
de a 
Using (A.1) along with the relationships 
a det (F*) tr(Vu) 
oa =tr 
de 20 
d E T 
—cof(F*) = tr(Vu)I- Vu 
de e=0 


a straightforward calculation yields 


o = u (Vu + Vu?) + Ar(Vu)I — (0 — 69)(3A + 20)I 


1 This Appendix is based on [1] 
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Introducing the infinitesimal strain tensor € = z(Vu + Vu?) the linearized stress tensor 
can be rewritten as 


c = 2pe + Atr(e)I — 3K (0 — 09)1 (A.2) 


Expression (A.2) for the stress tensor coincides with the generalized Duhamel-Neumann 
form of Hookes's law for isotropic behavior in the framework of linear thermoelasticity (cf. 
[159, Sec. 6.2]). Accordingly, A and u can be identified as the Lamé constants, K = A+ inu 
is the bulk modulus, and £ is the linear coefficient of thermal expansion. 


A.2. Lyapunov function 


Following [146] we consider the functional 
Li = € — 095! (A.3) 


where 0o is the reference temperature. Note that the total energy, €’, and the total entropy, 
S', have been defined in (4.39) and (4.41), respectively. The time derivative of (A.3) is 
given by 


d£ det, ds’ 
= ? dt 


Taking into account the results of Section 4.3.2, the last equation can be recast in the 
form 
dE. 


E =~ RESET) 655] - [8,87 


sau) 


1 
zN Qaa) 


boun 


=} p p.PNdA- N-QdA- & (Is. - 
OB OB 


= Fext — 0o ([S’, S)] + G.) 


0B 


Here, Pext = Jog p |p-tdA is the power expended on B by external Piola tractions 
t = PN acting on OB. Moreover, 


0 — 0o 
e = s dA A.4 
N (A4) 


is the entropy produced at the interface between the body and the environment. The 
environment is called thermally perfect in the sense of [146], if C, — 0. In the numerical 
examples we consider the case that Ge = 0 and P4 = 0, so that 


d£ 
dt 


= —-6o|S', S'] <0 


Then £’ is non-increasing over time and thus a natural Lyapunov function. 
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A.3. Convergence criteria 


To terminate the iterative solution procedure we use the norm of the residual vector as 
convergence criterion for both the (ME), scheme and the (M)g scheme. Since the (EM). 
scheme is capable to satisfy the balance of energy, we define an alternative termination 
criterion based on the discrete balance of energy: 


(Et — EB) - At (Phot + QM) | <e (A.5) 


Here, £ is the prescribed tolerance, E} is the discrete version of the total energy (4.39) at 
; —h a 
time tn, phex = Jon p Phys ‘thy dA is the power expended on the body by external 
nodal forces due to external tractions, and Qhint = — J. dB qi i dA is the heat flux across 
2 


the boundary. 
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B. Appendix to Chapter 5! 


B.1. Notes on the implementation 


We provide a short outline of the implementation of the present EME schemes. To this 


end, the residual vector associated with node a can be written in the form 


T 
R° -(nz, RO 


aT mn Ra, Re) 


(B.1) 


The nodal contributions to (B.1) emanate from weak form (5.65) by inserting the finite 


element approximations for the test functions leading to the expressions 


a) v 


RE 
2 
Rê = (m a= a Mn) +F} raSi VN") dV — 
B 


h — nb 
RZ = N° C Pr a 


3s B 


where, according to (5.66), 
Sr = 2 (Dcu^ — Oh, Don") 


1 
Qe zi (05, Kaigo V (s- )-- = Kigo V Oigo 
algo 


oh _ TI,(D-u") | N*(D;u)a 
er II; (D-n) 5 N®(D-n) 


! This Appendix is based on [2] 


(B.2) 
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Note that eventually (B.2) leads to a nonlinear algebraic system of equations, R^ — 0, 
a =1,...,N, for the determination of the nodal state variables qa, ||, va, ,,, and Ta, 4, 
along with the nodal quantities (D-u)., (D;7)q related to the projections II; (Du?) and 
IT, (D^), cf. (5.72). As usual, we apply Newton's method for that purpose. The consistent 
linearization of R^ can be performed in a straightforward way. In this connection, the 
projection formulas II; (D.-u^) and II„(D-n") do not pose any additional difficulties. For 
example, the linearization of the algorithmic temperature assumes the form 


A(D.u), ^ A(Dzm)a ) 
II; (Dzuh) II; (Dzh) 


AG ce NT ( 


Moreover, the material gradient of the algorithmic temperature is given by 


h | oh a (Dru)a (Dz7)a 
Vel, = elu" (Ira = mon 


leading to the corresponding linearized form 


D; Ja A(D-u)b 
A(V@ligs) = Oh, VN (a2 - Pe e 
(V Oigo) Eis V (st II;,(D,7/") II; (Dzu^) 
D. Ja (Dru)a A(D-n) 
_ oh N“ b 2 ( n N? 4 b 
Olgo V c II, (Dh) II (Dzu^) II; (D^) 


where 6° is the Kronecker delta. Concerning the termination criterion for the Newton 
iterations we apply the nonstandard form 


ED ael quA (/ 
0, 


which is made possible due to the energy consistency of the present (EME), schemes. In 
(B.3), € denotes the prescribed numerical tolerance for the Newton iterations. The left- 
hand side of (B.3) can be directly linked to the discrete balance of energy (5.80). Note that 
here, the standard contribution of the external tractions has been taken into account. 


<e (B.3) 


= h FE 
Em figs hd ney) 


q 
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C. Appendix to Chapter 6! 


C.1. Rotational symmetry 


Due to invariance under superposed rigid rotations the present problem satisfies symmetry 
properties (6.93). To see this, we consider a rotation tensor Q^ = exp(e£), expressed 
through the Rodrigues formula [24]. Here, € is a scalar parameter and 3 is a skew- 
symmetric tensor with associated vector £ € IR?, such that £a = € x a for all a € R?. 


Now, a superposed rigid motion gives rise to rotated nodal position vectors 
q; -Q'q,. (C.1) 
d NE 


Note that q? = q, and 3c me = €. Nodal pattern (C.1) of the rigidly rotated discrete 
system gives rise to the corresponding expression for the discrete deformation gradient 


(cf. (6.53)1) 
ph5—qpevN*-Q'Eh. (C.2) 
Furthermore, the corresponding right Cauchy-Green deformation tensor follows from 
(6.3) and is given by 
Che = (FEFE = ch, (C.3) 


where property (Q*)? = (Q*)-! of the rotation tensor has been taken into account. A 
frame-indifferent formulation implies that the internal energy density takes the form 
u(F, Tr, cz") = u(C,r, cz"). Thus, taking the derivative of u(F'*, r, cz") = 
u(C^, r, C,!) with respect to parameter £ and subsequently setting £ = 0 yields 


d h, h —1 
= |. u (Ein (Cpe) 


= Ügug : q; 9VN, 


d 

A 2=0 (C.4) 
= Opty : (£a, Q VN;) 

= (€ x qa) Š OpugV Nz 

— £: (qa x gu, VN?) , 


1 This Appendix is based on [4] 
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where u, is given by (6.52). Due to the arbitrariness of € € IR? (and vanishing body 
forces), the last equation coincides with symmetry condition (6.93)1, that is 


qa X Og, —0. 


Note that this condition holds as well for the specific choice 7 = u, due to (6.97). 


Similarly, symmetry condition (6.93)? results from the frame-indifference of the entropy 
density function 7(F, 7, C,') = (C, r, C or from (6.100); in the case T = 77. 
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